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FOREWORD 


NASTRAN® (NASA STRUCTURAL ANALYSIS) is a large, comprehensive, 
nonproprietary, general purpose finite element computer code for structural 
analysis which was developed under NASA sponsorship and became available to 
the public in late 1970. It can be obtained through COSMIC® (Computer 
Software Management and Information Center), Athens, Georgia, and is widely 
used by NASA, other government agencies, and industry. 

NASA currently provides continuing maintenance of NASTRAN through COSMIC. 
Because of the widespread interest in NASTRAN, and finite element methods in 
general, the Eighteenth NASTRAN Users' Colloquium was organized and held at 
The Red Lion Portland Center, Portland, Oregon on April 23-27, 1990. (Papers 
from previous colloquia held in 1971, 1972, 1973, 1975, 1976, 1977, 1978, 
1979, 1980, 1982, 1983, 1984, 1985, 1986, 1987, 1988 and 1989 are published in 
NASA Technical Memorandums X-2378, X-2637, X-2893, X-3278, X-3428, and NASA 
Conference Publications 2018, 2062, 2131, 2151, 2249, 2284, 2328, 2373, 2419, 
2481, 2505 and 3029.) The Eighteenth Colloquium .provides some comprehensive 
general papers on the application of finite element methods in engineering, 
comparisons with other approaches, unique applications, pre- and 
post-processing or auxiliary programs, and new methods of analysis with 
NASTRAN. 

Individuals actively engaged in the use of "inite elements or NASTRAN 
were invited to prepare papers for presentation at the Colloquium. These 
papers are included in this volume. No editorial review was provided by NASA 
or COSMIC; however, detailed instructions were provided each author to achieve 
reasonably consistent paper format and content. The opinions and data 
presented are the sole responsibility of the authors and their respective 
organizations. 


NASTRAN® and COSMIC© are registered trademarks of the National Aeronautics and 
Space Administration. 


i i i 

PRECEDING PAGE BLANK NOT FILMED 


' t 

RftSE , / / INTENTIONALLY BLANK 



CONTENTS 


Page 

FOREWORD i i i 


1. A NASTRAN TRAINER FOR DYNAMICS L 

by H. R. Grooms, P.J. Hinz and G. L. Commerford '/ 

(Rockwell International) 

2. PCI - A PATRAN-NASTRAN MODEL TRANSLATOR 14,' 

by T. J. Sheerer ~‘ 2 - 

(Chrysler Technologies Airborne Systems, Inc.) 

3. A GENERIC INTERFACE BETWEEN COSMIC/NASTRAN AND PATRAN 2 O 5 - 

by Paul N. Roschke, Prakit Premthamkorn , and James C. Maxwell 

(Texas A&M University) 

4. ACCURACY OF THE QUAD4 THICK SHELL ELEMENT 30. 

by William R. Case, Tiffany D. Bowles, Alicia K. Croft and 

Mark A. McGinnis 

(NASA/Goddard Space Flight Center) 


5. EIGENVALUE COMPUTATIONS WITH THE QUAD4 CONSISTENT-MASS MATRIX . . 61,- 

By Thomas A. Butler 

(Los Alamos National Laboratory) 


6. NASTRAN MIGRATION TO UNIX 1Z,- 

by Gordon C. Chan and Horace Q. Turner ‘ 

(UNISYS Corporation) 

7. OBTAINING EIGENSOLUTIONS FOR MULTIPLE FREQUENCY RANGES IN A 

SINGLE NASTRAN EXECUTION 77:,, 

by P. R. Pamidi and W. K. Brown 
(RPK Corporation) 


8. RANDOM VIBRATION ANALYSIS OF SPACE FLIGHT HARDWARE USING NASTRAN 90 
by S. K. Thampi and S. N. Vidyasagar 
(GE Government Services) 


9. OBTAINING AN EQUIVALENT BEAM 107; 

by Thomas G. Butler 
(Butler Analyses) 

10. LOW VELOCITY IMPACT ANALYSIS WITH NASTRAN 115_ 

by Daniel A. Trowbridge 
(Anal ex Corporation) 

and Joseph E. Grady and Robert A. Aiello 
(NASA Lewis Research Center) 


v 


PRECEDING PAGE BLANK NOT FILMED 


mm 



CONTENTS 

(Continued) 


Page 


11. TRANSITIONING OF POWER FLOW IN BEAM MODELS WITH BENDS 135 

by Stephen A. Hambric yy 

(David Taylor Research Center) 

12. COUPLED BE/FE/BE APPROACH FOR SCATTERING FROM FLUID-FILLED 

STRUCTURES 150y 

by Gordon C. Everstine and Raymond S. Cheng 'h 

(David Taylor Research Center) 

13. MONITORING OF RITZ MODAL GENERATION 165 

by Mladen Chargin 

(NASA Ames Research Center) 
and Thomas G. Butler 

(Butler Analyses) 


' £ '«•*• ■* 





A NASTRAN TRAINER FOR DYNAMICS 

H.R. Grooms, P.J. Hinz, and G.L. Commerford 
Rockwell International 


N 9 0 - 24^ /8 8 


OVERVIEW 

As the use of the finite element method proliferates , the need for training 
becomes more and more pronounced. An automated tool to familiarize engineers with 
static solutions has been developed and used. This tool (Figure 1) is part of an 
overall structural analysis/expert training system (ref. 1). Experiences with this 
tool and comments from users (ref. 2) have underlined the need for a dynamic version 

of the trainer. 

{ r •• 

This paper presents an automated training tool that engineers can use to master 
the application of NASTRAN to dynamic problems. The paper consists of the following 
sections: 


• Overview 

• Background 

• Existing Programs 

• Scope, Purpose, and System Organization 

• Example Problems 

• Conclusions 

• References 

Example problems have been selected to make classical solutions available for 
comparison. These comparisons can be used to evaluate the solution. 


BACKGROUND 

The solution of dynamic problems involves some complications that do not exist 
with static problems: 

• How many degrees of freedom should be retained for the eigenvalue solution? 

• Which discrete mass items are so large or important that they should be 
retained for eigenvalue solution? 

• How many frequencies and mode shapes are needed and to what accuracy? 

An engineer may think that most of the mass associated with a structure can be 
traced to the structural members themselves; this is not necessarily true. With many 
aircraft and spacecraft, the nonstructural masses (e.g., hydraulic lines, fuel tanks, 
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environmental control equipment, etc.) have a pronounced influence on the overall 
mass distribution and may have the greatest dynamic effect. 

The example problems have distributed masses and lumped masses that the user 
must consider in the solution approach. These examples help the user develop judgment 
when deciding on the number and the particular degrees of freedom to be retained, and 
on how to discretize the distributed mass. 


EXISTING PROGRAMS 

Various researchers have developed computer programs for structural analysis and 
design applications. Ginsburg (ref. 3) addressed computer literacy, while Woodward 
and Morris discussed improved productivity through interactive processing (ref. 4). 
Wilson and Holt (ref. 5) developed a system for computer-assisted learning in 
structural engineering. Sadd and Rolph (ref. 6) described the various ways in which 
design engineers could be trained to use the finite element method. Self-adapting 
menus for CAD software are covered by Ginsburg (ref. 7). 

Bykat (ref. 8) is developing a system that will have features for training, 
analysis control, and interrogation. 


SCOPE, PURPOSE, AND SYSTEM ORGANIZATION 

The NASTRAN trainer was designed to be a stand-alone tool. The trainer is user 
friendly — a knowledge of job control language or the operating system is not 
required. A user can sit down at a terminal and, in very little time, start solving 
an example problem. The trainer is organized so that a user must complete the static 
problems (ref. 2) before the dynamics problems can be accessed. This organization 
prevents a user who has no familiarity with the finite element method from starting 
with the dynamics section. 

The trainer is organized into three main modules: (1) overview, (2) user’s 
guide, and (3) problem set. Figure 2 shows some details of each module. The user 
accesses these modules by using the primary menu. More details of the NASTRAN 
environment sections are given in Figure 3. 


EXAMPLE PROBLEMS 

The example problems, shown in Figures 4 through 11 and summarized in Table 1, 
become progressively more difficult to solve. The first problem is a simply supported 
beam with a single lumped mass at the center. 

There are various courses and classes to instruct engineers in solving dynamics 
problems. These courses usually emphasize the theory. A vital part of solving any 
large dynamics problems is deciding how many and which degrees of freedom should be 
retained for the eigenvalue solution. This is usually a matter of judgment, and it 
takes solving many problems to develop this judgment. 
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Example 2 was solved using three different approaches. The user was trying to 
answer some fundamental questions that must be addressed every time a dynamics 
problem is solved using the finite element method: 

• Is the model fine enough? 

• Have the distributed masses been lumped into enough locations? 

• Have enough degrees of freedom been retained in the eigensolution? 

Figure 12 summarizes the different approaches. Table 2 compares the computed 
three lowest natural frequencies with the exact results. 


CONCLUSIONS 

An automated training tool that helps engineers become familiar with using 
NASTRAN to solve dynamic problems has been presented. The tool allows the user to 
proceed at his own pace by using a set of eight example problems. The examples were 
selected so that classical solutions are available and displayed, enabling the user 
to make comparisons. 
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Table I. Example Problems 


Example 

Description 

Significant Features 


Beam simply supported on both ends with lumped mass in middle 

Motion in one plane only, lumped mass only 


Beam simply supported on both ends with uniformly distributed 
mass 

Motion in one plane only, distributed mass 


Beam fixed on one end with a lumped mass at the free end 

Motion in any direction, lumped mass only 


Beam fixed on one end with a uniformly distributed mass 

Motion in any direction with uniformly distributed mass 


Rectangular plate clamped on one edge, all other edges free with 
a uniformly distributed mass 

Plate bending with distributed mass 


Rectangular plate, free-free with uniformly distributed mass 

Free-free (implies six modes with zero frequency) 

7 

Two beams connected by springs, each with distributed and 
lumped mass 

Multibody problem, free-free 

8 

Problem 7 with a forcing function added 

Forcing function 


Table 2. Comparison of Natural Frequencies for Example 2 


Mode 

Exact 

Solution 

First Approach 

Second 

Third 

i 

9.870 

9.867 

9.869 

9.872 

2 

39.48 

39.19 

39.47 

39.74 

3 

88.83 

83.21 

88.66 

93.62 
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Figure 2. Organization of NASTRAN Trainer 
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Organization of Program 
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Three Approaches to Example 2 
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PCI - A PATRAN-NASTRAN MODEL TRANSLATOR 



T.J. SHEERER 

CHRYSLER TECHNOLOGIES AIRBORNE SYSTEMS, 


INC. 


WACO, TX 


INTRODUCTION 


The existence of several derivative versions of NASTRAN, which 
differ significantly in element definitions and result formulation, 
has caused some difficulties in the interface between NASTRAN and 
pre- and post- processors such as PATRAN or SUPERTAB. In 
particular, the PATRAN-COSMIC/NASTRAN interface provided by PDA 
Engineering has not been updated at the same rate as the equivalent 
interface with MSC/NASTRAN, and has significantly less 
capabilities. Many model entities supported by both PATRAN and 
COSMIC /NASTRAN are not supported by the translator. The well- 
documented PATRAN neutral file, which is now supported by several 
other vendors, has provided a means for the user to create his own 
interface program for model translation, while it has also been 
possible to pass results form NASTRAN to PATRAN with a user-written 
program using 0UTPUT2 statements and format information from the 
Patran Users' Manual. In recent years PDA engineering has provided 
as part of PATRAN a library of subroutines known as gateway 
utilities, which extract data directly from the PATRAN database 
file and which can be called from a FORTRAN program. As this 
eliminates the task of reading the neutral file the work of 
creating a translator can be produced by the user with little 
effort. This has been done with the object of producing a PATRAN- 
COSMIC NASTRAN model translator comparable in scope to the PATRAN- 
MSC/NASTRAN translator, but also allowing a greater degree of user 
control than is found therein. The different parts of the program 
were developed in several locations, as the counterpart to a 
results translator developed for Texas Instruments by Texas A and 
M University, which also emphasizes flexibility. Both these 
programs are public domain programs under the terms of the 
development agreement between Texas Instruments and Texas A.M.U., 
and enhancements developed at Chrysler have also been passed to 
T.A.M.U. 
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PROGRAM STRUCTURE 


PCI supports a range of elements comparable with PATNAS and 
significantly greater that PATCOS. The structure of the program 
is such that, in effect, it supports elements not currently extant. 
In the PATCOS program element types are hard-wired, and if a 
different NASTRAN card is required a text editor must be used to 
change the bulk data file. If, for example, the QUAD4 element is 
required in a model, all elements of this type must be edited from 
type QUAD2 to QUAD4, since PATCOS does not currently support the 
QUAD4. Similar difficulties arise with other elements. The PATNAS 
translator does support a wide range of elements identified by 
configuration codes, but the equivalence between elements and 
configuration code is controlled from within PATRAN, and the user 
can not translate to an element developed in-house, or a newly 
introduced element, or simply an element not supported by the 
translator. PCI uses a different approach by having a user text 
file of twenty NASTRAN element names for each element 
configuration. This file may be edited by the user to assign 
whatever correspondence he desires between configuration code and 
NASTRAN element type. Since most element cards in NASTRAN follow 
the same pattern of (NAME, PID/MID, NODES) it is not generally 
necessary to write a new subroutine when adding an element. 
Frequently the only parameter changed is field one of the NASTRAN 
data card, and this is controlled from a text file. Because 
NASTRAN and PATRAN use a different numbering sequence for higher 
elements it is convenient to use several subroutines to write 
NASTRAN elements, but frequently only one routine is required for 
each shape/node combination. All 3-node shells, for example, are 
written by the same routine. The PATRAN configuration field is 
used to select an element name from a range of twenty for each 
configuration of shape/nodes. The user-generated text file 
ETYPES.DAT contains these names for configurations 1 through 20. 
Any elements having the default configuration of zero in the PATRAN 
database are assigned the value 1 in the translator. 

Table (1) shows the basic PATRAN element types supported by 
PCI and the NASTRAN elements obtainable form them. For comparison 
the elements supported by PATCOS and PATNAS are also shown, note 
that, in general, failure by PCI to support elements listed for 
PATNAS is because they are not supported in COSMIC /NASTRAN. In 
practice, alteration of the text file of element configurations 
will allow support of these elements if and when they become 
available 
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TABLE 1 : GEOMETRIC ELEMENT TRANSLATION 




CAPABILITIES 

OF PATCOS, PATNAS AND PCJ 



shape/nodes 

PATCOS 

PATNAS 

PCI 


BAR/ 2 

CBAR 

CBAR 

CBEAM 

CBAR 




CROD 

CROD 


TRIA/3 

CTRIA2 

CTRIA1 

CTRIAl 




CTRIA2 

CTRIA2 




CTRIA3 

CTRBSC 




CTRIARG 

CTRIARG 

CTRMEM 

CTRPLT 


QUAD/ 4 

CQUAD2 

CQUAD1 

CQUAD1 




CQUAD2 

CQUAD2 




CQUAD4 

CQUAD4 




CQDMEM1 

CQDMEM1 




CQDMEM2 

CQDMEM2 




CSHEAR 

CSHEAR 

CTWIST 

CQDPLT 

CQDMEM 




CTRAPARG 

CTRAPAX 


TRIA/6 

CTRIM6 

CTRIA6 

CTRIM6 




CTRIAX6 

CTRSHL 

CTRPLT1 


QUAD/ 8 

CIS2D8 

CQUAD8 

CIS2D8 


HEX/ 8 

CIHEX1 

CHEXA 

CIHEX1 




CHEX8 

CHEXA1 

CHEXA2 


QUAD/ 20 

CIHEX2 

CHEXA 

CIHEX1 




CHEX8 

CHEXA1 

CHEXA2 



QUAD/ 3 2 


CIHEX3 
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PCI also supports the SCALAR and damping elements CELAS2 and 
CDAMP2 by generating a scalar element for each degree of freedom 
specified in a PATRAN SPRING element. The concentrated mass 
element C0NM2 is obtainable from the PATRAN MASS directive. MPCs 
and rigid elements are supported. Table (2) summarizes the support 
for these elements. Node translation with embedded SPCs is fully 
supported . 


TABLE 2: 

SCALAR, DAMPER AND MASS 

ELEMENT SUPPORT 

PATRAN DIRECTIVE 

NATRAN 

CARD WRITTEN BY PCI 


MPC 


MPC 

MPC , RROD 


CRIGIDR 

MPC , RBAR 


CRIGD1 

MPC, RBE1 


CRBE1 

MPC , RBE2 


CRBE2 

MPC, RBE3 


CRBE3 

BAR/2/n OR 

SPRING 

CELAS2 

BAR/ 2 /n OR 

MASS 

C0NM2 

BAR/2/n 


CDAMP2 


LOAD AND CONSTRAINT TRANSLATION: 

PCI supports constraint (SPC), nodal force and pressure 
translation. FORCE and SPC cards are translated in an element- 
dependent manner as shown in Table ( 3 ) . Only normal element 
pressures are supported. PCI will select the appropriate PLOAD 
card for the element type. In many circumstances involving higher 
order elements several PLOAD cards must be generated for a single 
PATRAN pressure load . 
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TABLE 3: 

PCI SUPPORT OF LOADS AND SPCS 


NASTRAN CARD 

USED IN ELEMENTS 

SUPPORT 

PL0AD4 

CQUAD4 

YES 

PL0AD2 

OTHER 1ST ORDER SHELL 

YES 

PL0AD3 

ISOPARAMETRIC SOLID 

NOT YET 

PLOAD 

2ND ORDER SHELL 

YES 


COORDINATE SYSTEMS: 


The various PATRAN coordinate systems are translated to C0RD2 
cards in NASTRAN, as in the PATNAS translator. Additionally, 
PATRAN nodes, which are stored in the database in the basic 
coordinate system, have their locations output to NASTRAN in the 
local system associated with the nodes. This is of great 
importance if, for example, constraints are to be applied in a 
local system. 


The PATRAN database includes, associated with each coordinate 
system, a 3x3 matrix T such that 





where the suffixes 1, b, denote local and basic coordinate systems 
respectively. The translator inverts the matrix to produce the 
matrix T' 1 where 



and prints the local values to the NASTRAN data deck. 


PROPERTY GENERATION: 

It is generally no more difficult to enter property and 
material cards in NASTRAN that in PATRAN. For this reason property 
generation has not been implemented in PCI . 
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CONCLUSIONS: 


The amount of programming required to develop a PATRAN-NASTRAN 
translator was surprisingly small. The approach taken produced a 
highly flexible translator comparable with the PATNAS translator 
and superior to the PATCOS translator. The coding required varied 
from around ten lines for a shell element to around thirty for a 
bar element, and the time required to add a feature to the program 
is typically less than an hour . The use of a lookup table for 
element names makes the translator also applicable to other 
versions of NASTRAN. The saving in time as a result of using PDA s 
Gateway utilities was considerable. 

During the writing of the program it became apparent that, 
with a somewhat more complex structure, it would be possible to 
extend the element data file to contain all data required to define 
the translation from PATRAN to NASTRAN by mapping of data between 
formats. Similar data files on property, material and grid formats 
would produce a completely universal translator from PATRAN to any 
FEA program, or indeed any CAE system. 
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1 ^ a GENERIC INTERFACE BETWEEN COSMIC/NASTRAN AND PATRAN 

Paul N. Roschke, Prakit Premthamkorn, and James C. Maxwell 
Texas A&M University 

SUMMARY 

Despite its powerful analytical capabilities, COSMIC/NASTRAN lacks adequate 
post-processing adroitness. PATRAN 1 , on the other hand is widely accepted for its 
graphical capabilities. A nonproprietary, public domain code mnemonically titled CPI (for 
COSMIC/NASTRAN-PATRAN Interface) is designed to manipulate a large number of 
files rapidly and efficiently between the two parent codes. In addition to PATRAN's results 
file preparation, CPI also prepares PATRAN's P/PLOT data files for xy plotting. The user 
is prompted for necessary information during an interactive session. Current 
implementation supports NASTRAN’s displacement approach including the following rigid 
formats: (1) static analysis, (2) normal modal analysis, (3) direct transient response, and 
(4) modal transient response. A wide variety of data blocks are also supported. Error 
trapping is given special consideration. A sample session with CPI illustrates its simplicity 
and ease of use. 


INTRODUCTION 

Overview 

The standard gateway that interfaces COSMIC/NASTRAN's analysis results to 
PATRAN's post-processing makes use of NASTRAN's FORTRAN-written results files. 
These files can be requested via DMAP ALTER's OUTPUT2 statement in NASTRAN's 
executive data deck. They contain data in mixed ASCII and binary code format. However 
they can not be used as direct input to PATRAN. Similarly, PATRAN also supports 
communications with external codes via specially formatted results files. Format of these 
files is predetermined according to PATRAN and differs for each data type. Generally, 
they can be categorized into three groups according to their formats: (1) nodal results files’ 
(2) element results files, and (3) beam results files. The number of results files can be as 
many as required. Therefore, in order to interface COSMIC/NASTRAN to PATRAN for 
post-processing purposes, an interface that is capable of translating from NASTRAN's 
analysis results to PATRAN-recognizable files is required. 

Fig. 1 shows code and file relationships among the analysis and post-processing 
modules. The flow begins with the input of NASTRAN's analysis data deck. DMAP 
ALTER statements must be included in the executive control deck in order to obtain a 
FORTRAN-written results file. Generally, this file is assigned an extension of .PAT. This 
NASTRAN results file is then translated by the COSMIC/NASTRAN-PATRAN Interface 
(CPI). Output from CPI is a group of PATRAN-compatible results files and/or a 
P/PLOT-compatible data file. These files are ready for PATRAN post-processing. 


l 


PATRAN is a trademark of PDA Engineering 
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PATRAN also requires access to a neutral file containing geometrical properties of the 
model. This file can be obtained via the COSPAT 2 translator, but not by means of CPI. 



Fig. 1 Code and File Relationships 


Program Description 

COSMIC/NASTRAN-PATRAN Interface (CPI) translates COS MIC/NASTR AN s 
FORTRAN-written results files to PATRAN compatible results files. These results files 
can be requested via the OUTPUT2 instruction in the executive control deck. CPI 
provides the operation in two modes: 

Mode 1 - Whole Model Translation: CPI translates all data blocks contained in the 

NASTRAN results file (*.PAT) and creates PATRAN compatible results files 
corresponding to data blocks found in the input data file. The user is prompted for the 
prefix name of an output file. Different prefixes allow the user to distinguish between 
groups of output files when many results files are translated in a single execution Up to 6 
characters per prefix is acceptable. Created output files and translated data blocks are 
summarized on the screen during CPI's execution. 

Mode 2 - XY-Plot Data Preparation: CPI creates data for the specified node or element 
from a set of data blocks. CPI prompts the user for necessaiy information, i.e node 
number and component number for nodal data, and searches the input file for any data on 


2 COSPAT is a COSMIC /NASTRAN to PATRAN interface, developed by COSMIC, University of Georgia, 
Athens, Georgia. 
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this node. Each data set is written to PATRAN's P/PLOT compatible data file. CPI 
names this file XYPLOT.DAT. When this command is selected, CPI automatically reviews 
the data blocks available in the current input file. The user is then prompted to enter the 
number corresponding to the data block name shown above this prompt. 

Currently, CPI supports NASTRAN's displacement approach including, but not 
limited to, the following rigid formats: 

(1) Static Analysis (Rigid Format 1) 

(21 Normal Modal Analysis (Rigid Format 3) 

(3) Direct Transient Response (Rigid Format 9) 

(4) Modal Transient Response (Rigid Format 12) 

Table 1 shows data blocks supported by CPI. It should be noted that CPI supports 
any rigid format as long as the data blocks listed are encountered. Rigid formats named 
above are only a guide to the rigid formats most frequently giving rise to these data blocks. 
CPI recognizes data blocks, not rigid formats. Table 2 shows elements currently supported 
by CPI. J 


Table 1. Data Blocks Supported by CPI 


Data Block 

Content 

OUGV1 

Nodal Displacements 

HOUGVl 

Nodal Temperature 

OQG1 

Single Point Constraint 

OPG1 

Load Vectors 

OPHIG 

Eigenvectors 

OUPV1 

Displacement Vectors 

ONRGY1 

Strain Energy 

OESl(X) 

Element Stress 

OEF1 

Element Forces 

OGPFB1 

Grid Point Balance Forces 


Table 2. Elements Supported by CPI 


Element Name 

ID Number 

Element Name 

ID Number 

CROD 

1 

CQDPLT 

15 

CQDMEM2 

63 

CBEAM 

2 

CTRIA2 

17 

CTUBE 

3 

CQUAD2 

18 

CTRIA1 

6 

CQUAD1 

19 

CTRBSC 

7 

CBAR 

34 

CTRPLT 

8 

CTRIARG 

36 

CTRMEM 

9 

CTRAPRG 

37 

CONROD 

10 

CQDMEM1 

62 

CTETRA 

39 

CWEDGE 

40 

CHEXA1 

41 

CHEXA2 

42 

CIHEX1 

65 

CIHEX2 

66 

CIHEX3 

67 
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PROGRAM USAGE AND SAMPLE SESSION 

COSMIC/NASTRAN-PATRAN INTERFACE (CPI) is an interactive program. A 
user can manipulate the program by selecting from menu commands and answering 
questions. No special commands are needed. CPI also provides extensive error trapping to 
ensure appropriate input and output. 

The following sample session (adapted from Ref. 5) illustrates execution of CPI for 
transient dynamic analysis results in which a series of time steps is involved. A simple three 
node model which is 350 cm. high, and has a cross sectional area of 6.45 cm. 2 (see Fig. 2), is 
used to model a rocket trajectory. Time versus loading of the rocket is shown in Fig. 3. 


node 



Fig. 2 Sample 2 - Geometric Properties 


Load (N) 


| \ Time (sec) 

1.1 

1.125 

Fig. 3 Sample 2 - Load vs. Time 

After invocation, CPI's opening header appears on the screen and CPI prompts the 
user for the name of the NASTRAN results file to be translated. The user then enters the 
NASTRAN results file name with its extension. Alternatively, the user can leave the 
translator and return to the operating system by entering "EXIT' (or "exit" or "E" or simply 
"e"). Here, the NASTRAN results file name is assigned the name "ROCKET.PAT." 



0.0875 
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****************************************************** 

** 


** 

** 

COSMIC/NASTRAN - PATRAN INTERFACE 

** 

** 


** 

* * 

** CPI** 

** 

* * 


** 

* * 

TEXAS A&M UNIVERSITY 

** 

** 

& 

** 

** 

TEXAS INSTRUMENTS 

** 

** 


* * 

** 

(JUNE 1989) 

* * 

** 


** 

********* 


******* 

ENTER NASTRAN RESULTS FILE NAME OR "EXIT" 
>ROCKET . PAT 



Once the results file has been named, CPI reads the results file header and echoes 
it to the screen for verification. Next a menu and a prompt for selection appear on the 
screen. The response to CPI should be one of the following: 

> 1 Select 1 to translate the entire file and create appropriate PATRAN-compatible 
results files. 

> 2 Select 2 to translate the data and create PATRAN's P/PLOT-compatible data files 
which contain requested nodal or element data. 

> 3 Select 3 to specify a new NASTRAN results file. 

> 4 Select 4 to exit CPI and return to the operating system. 

If the first option is chosen, translation proceeds as follows: 


DATE : 7/27/89 

HEADER: NASTRAN FORT TAPE ID CODE - 
LABEL : XXXXXXXX 


>>>>>>> ENTER YOUR SELECTION 


1) WHOLE MODEL TRANSLATOR 

2) XY-PLOT DATA 

3) ENTER NEW RESULTS FILE 

4) EXIT 

> 1 

ENTER FILENAME PREFIX (UP TO 6 CHARACTERS) OR [RETURN] 
(EXAMPLE DEFAULT FORMAT: S I .DIS) 

> ROCKET 

IS "ROCKET" THE CORRECT PREFIX? [Y] 

><RETURN> 
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At this point the entire model is translated with the creation of the 51 files shown on 
the screen during execution. Each file corresponds to a new time step and a new geometric 
location of the rocket. It is apparent that running each file with PATRAN to obtain 
specific data for a given node and component can become tedious. Therefore, a routine 
has been incorporated into CPI to allow the user to create a file that contains only user- 
specified data available for plotting with P/PLOT. 

When the second option is chosen, CPI reviews data block names existing in the 
current results file. The user is asked to select a data block name. CPI classifies data 
blocks into two groups: (1) nodal data, and (2) element data. If the selected data block 
contains nodal data, CPI prompts for a node number and a component number (1-6). In 

general, the six components of nodal data are: 

Component 1: X-direction translation vector 
Component 2: Y-direction translation vector 
Component 3: Z-direction translation vector 
Component 4: X-direction rotational vector 
Component 5: Y-direction rotational vector 
Component 6: Z-direction rotational vector 

If the selected data block is an element data block, the user is prompted for an element 
number and a column number. 

Next, CPI searches for the requested data and, upon completion of gathering the 
data, requests a data header extension of the current XY-plot data. Up to 38 characters is 
acceptable. Execution proceeds as follows: 
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»»»> ENTER YOUR SELECTION «««< 


1) WHOLE MODEL TRANSLATOR 

2) XY-PLOT DATA 

3) ENTER NEW RESULTS FILE 

4) EXIT 


DATA BLOCK REVIEW 


1) OUPV1 


> 1 
> 2 
> 2 


2 - 


ENTER DESIRED BLOCKNAME NUMBER 
(OR ''O'* TO EXIT) 

ENTER NODE NUMBER 

WHICH COMPONENT (1-6)? 

ENTER PLOT TITLE EXTENSION FOR DATA SET: 
DEFAULT TITLE: YDATA, DATA SET: 1; NODE: 


> Y-TRAJECTORY OF ROCKET 

%%%%% NUMBER OF DATA WRITTEN = 51 %%%%% 


1 

2 ; 


COMPONENT: 


DATA BLOCK REVIEW 


1) OUPV1 


ENTER DESIRED BLOCKNAME NUMBER 
(OR "0" TO EXIT) 

> 0 


»»»> ENTER YOUR SELECTION «<«« 


1) WHOLE MODEL TRANSLATOR 

2) XY-PLOT DATA 

3) ENTER NEW RESULTS FILE 

4) EXIT 

> 4 

CPI Execution Completed 


At this point a file called XYPLOT.DAT exists and upon input of this file into P/PLOT, 
data is readily graphed for the Y coordinates (component 2) of node 2 of the rocket for 
each of the 51 time steps (Fig. 4). 
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Fig. 4 Displacement of Node 2 - Component 2 


OUTPUT FILES 
Results Translation 

CPI creates several results files which provide a direct avenue between COS- 
MIC/NASTRAN's analysis algorithms and PATRAN’s post-processing capabilities. 
COSMIC/NASTRAN's input file (commonly called the input data deck) provides the 
name of a binary file which CPI interpolates. For example, EXAMPLE.NAS (the input 
data deck) yields EXAMPLE.PAT (the binary file to be translated to PATRAN- 
compatible results files). Creation of the binary file is accomplished by inclusion of the 
appropriate ALTER statement in the executive control deck of COSMIC/NASTRAN's 
input data deck. Only filename extensions are changed by execution of COS- 
MIC/NASTRAN. EXAMPLE.PAT is then translated by CPI into appropriate PATRAN- 
compatible results files. A description of each file follows. 

All CPI output files are given the name 'Xiilnnn.EXT' where X is either an 'S' or an 
'M' depending on whether the following two digit number is either a subcase or a mode 
number, respectively, and EXT is the file extension named with respect to data type. 
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XY-Plot 


Selecting "XY-PLOT DATA" at CPI's main menu initiates a prompt screen which 
lists the various data blocks found in COSMIC/NASTRAN's binary output file. CPI then 
invites the user to specify desired data blocks containing either nodal or element data that 
are to be written to XYPLOT.DAT, which is the input file for P/PLOT. This file can 
contain one or more data sets depending on how many times the user requests that CPI 
write to this file. CPI writes only y-data to XYPLOT.DAT since the user defines an initial 
x and delta-x when executing P/PLOT. A description of XYPLOT.DAT's format is given 
below. 

(1) Nodal Data Blocks 

The first line written for each data block contains the default title: 

•'YDATA, DATA SET: iii; NODE nnnn; COMPONENT j - » 

and may be appended to allow for more descriptive titles. This appendage plus the default 
title must not exceed 80 characters. If more than 80 characters are specified, those beyond 
column 80 are truncated. This allows the user a suffix of anywhere from 33 to 38 characters 
depending on the number of digits contained in the default title. 

Each line thereafter, until either a *Z (end of file) or another title is encountered, 
contains nodal data for components 1, 2, 3, 4, 5, or 6, i.e., X, Y, Z, © x , e y , or e z , 
respectively, as requested from CPI’s prompting. 

(2) Element Data Blocks 

As with nodal data blocks, the first line written for an element data block contains a 
default title and takes the form: 

'•YDATA, DATA SET: iii; ELEMENT nnnn; COLUMN j - " 

It may also be appended as described above, but again, the default title plus the appendage 
must not exceed 80 characters in length or truncation of characters beyond column 80 
occurs. As with the nodal data case, a suffix of 33 to 38 characters is allowed. 

The remaining lines associated with the current data set contain element data 
available for plotting. The data type is chosen by entry of a column number at CPI's 
prompt. 


SOFTWARE AND HARDWARE REQUIREMENTS 

CPI source code was written in standard FORTRAN 77; however, some special 
features of VAX/VMS FORTRAN were implemented. CPI will run on any VAX machine 
with a VMS operating system and a FORTRAN compiler. The total length of the source 
code is approximately 1200 lines, which requires 104 blocks of disk space. 


28 



REFERENCES 


(1) NASTRAN User's Manual, Computing Services Annex, University of Georgia, Athens, 
Georgia, 1986. 

(2) NASTRAN Programmer's Manual, Computing Services Annex, University of Georgia, 
Athens, Georgia, 1978. 

(3) PATRAN Plus User Manual, Release 2.3, PDA Engineering, Costa Mesa, California, 
1988. 

(4) Schaeffer, G. H., MSC/NASTRAN Primer, Schaeffer Analysis Inc., Mont Vernon, New 
Hampshire, 1979. 

(5) MSC/NASTRAN Demonstration Problem Manual, The Macneal-Schwendler 
Corporation, Los Angeles, California, 1983. 

(6) Nyhoff, L., Leestina, S., FORTRAN 77 for Engineers and Scientists, Macmillan 
Publishing Company, New York, 1988. 


29 



N9 0-246^1 

1 3 ACCURACY OF THE QUAD4 THICK SHELL ELEMENT 

/>) j\^ , William. R. Case, Tiffany. D. Bowles, Alicia K. Croft* 

and Mark. A. McGinnis 
NASA/Goddard Space Flight Center 


SUMMARY 


The accuracy of the relatively new QUAD4 thick shell element is assessed via 
comparison with a theoretical solution for thick homogeneous and honeycomb 
flat simply supported plates under the action of a uniform pressure load. The 
theoretical thick plate solution is based on the theory developed by Reissner 
and includes the effects of transverse shear flexibility which are not 
included in the thin plate solutions based on Kirchoff plate theory. In 
addition, the QUAD4 is assessed using a set of finite element test problems 
developed by the MacNeal-Schwendler Corp . (MSC) . Comparison of the COSMIC 
QUAD4 element as well as those from MSC and Universal Analytics, Inc. (UAI) 
for these test problems is presented. The current COSMIC QUAD4 element is 
shown to have excellent comparison with both the theoretical solutions and 
also those from the two commercial versions of NASTRAN that it was compared 
to . 


INTRODUCTION 


The QUAD4 thick shell element, added to NASTRAN in 1987, is one of the most 
important additions to the program since the original writing of the code. 

The deficiencies of the original QUADl and QUAD2 quadrilateral shell elements 
have been recognized for years and have been reported in the literature. At 
the Goddard Space Flight Center (GSFC) , the quadrilateral shell element is in 
use in virtaully all structural analyses of our spaceraft and related 
hardware. Typical applications are for the modelling of cylindrical shells 
and flat plates made of honeycomb or machined, lightweighted, metal that make 
up the structure of spacecraft and scientific instruments. In some cases 
these models require that the effects of transverse shear flexibility be 
included due to their thickness. The QUAD4 element includes these effects 
and, in addition, has an improved isoparametric membrane capability for 
in-plane loading. 

The purpose of the study reported herein is to assess the accuracy of the 
QUAD4 element in modelling a variety of situations involving both solid 
cross-section plates as well as those constructed of honeycomb. Three goals 
of the study were to determine: 

a) what is the rate of convergence to the theoretical solution as the 
mesh is refined; 


* Currently with Swales and Associates, Inc. 
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b) whether the element exhibits sensitivity to aspect ratios 
significantly different than 1.0; 

c) how the element behaves in a wide variety of modelling 
situations, such as those included in the MSC element test library 
(discussed below) . 

The first two questions were addressed in the same manner as several other 
studies reported by one of the authors in prior NASTRAN colloquia (references 
1 and 2). The procedure used in those studies, and followed here also, is to 
isolate the effects of mesh refinement and aspect ratio. That is, the mesh 
refinement study is done using elements with an aspect ratio of 1.0. Then, 
once a fine enough mesh has been reached such that the errors are small, the 
effects of aspect ratio can be investigated by keeping the mesh the same (i.e. 
same number of elements) and varying the overall dimensions of the problem, 
thus resulting in each element aspect ratio changing. Obviously, in order to 
accomplish this latter step there must be a theoretical solution (or some 
other equally acceptable comparison solution) to the problem with which to 
compare the finite element model results. This is needed since, at each step, 
a problem of different dimensions (and therefore different theoretical 
solution) is being modelled. 

The above tests are important in that they show the rate of convergence toward 
the theoretical solution as the mesh is refined. Those tests, however, are 
not sufficient to completely test the accuracy of a finite element since they 
do not test irregular geometries, or a variety of loadings or material 
properties. The MSC has developed a comprehensive set of problems for testing 
finite elements in a variety of situations (reference 3). The library of 
problems consists of 15 test problems for the QUAD4 element that cover all of 
the parameters mentioned above. A test of the COSMIC QUAD4 using these 
elements was reported at the 17th NASTRAN Users Colloquium in 1987 by Victoria 
Tischler of the Air Force Wright Aeronautical Laboratories (AFSC) at Wright 
Patterson Air Force Base, Ohio, but was not included in the formal 
proceedings. Due to the fact that it was not included in the formal 
proceedings, and also due to the fact that errors in the QUAD4 code for 
nonhomogene ous plates (to be discussed later) have been corrected, the results 
of our testing of the latest version of the element with the MSC library are 
include herein. 


RESULTS OF MESH AND ASPECT RATIO STUDY 


For the mesh and aspect ratio study a theoretical comparison solution is 
highly desirable. Since the effects of transverse shear flexibility are 
included in the QUAD4 element formulation, a theoretical solution for 
moderately thick plates, based on Reissner (or Mindlin) thick plate theory is 
also desirable. Such a solution is given in references 4 and 5 for 
rectangular simply supported thick plates under the action of a pressure load. 
Thus, this problem was used for the mesh and aspect ratio portions of the 
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study. Figure 1 defines the geometry, coordinate system, boundary conditions 
and loading for the rectangular plate. The thickness indicates a moderately 
thick plate of length to thickness ratio of .20. The effect of transverse 
shear flexibility is only approximately 1% on the maximum displacement but is 
important in discerning the quality of the convergence of the finite element 
results to the exact theoretical solution. By exact is meant the 
theoretical basis for the QUAD4 element, which is expressed in the Reissner 
thick plate theory. Figure 2 shows the finite element mesh geometry used in 
the mesh and aspect ratio studies. Due to symmetry only one quarter of the 
plate was modelled The 4x4 mesh shown on figure 2 is an example only; the 
mesh was varied during the mesh study. 

Figures 3a - 3c show characteristics of the theoretical solution. As indicated 
in figure 3a the central displacement solution is represented as an infinite 
series of hyperbolic functions. A FORTRAN computer program was written to 
compute the theoretical solutions for displacements (using the series shown) 
as well as stresses (solution not shown). As m gets large, where m is the 
number of terms included in the series, the hyperbolic functions tend to 
overflow the exponent range of the computer. This does not indicate a problem 
with the series shown, as the hyperbolic functions appear in both the 
numerator and denominator and their ratio is numerically stable. However, in 
separately evaluating the numerator and denominator the overflow problem was 
encountered. In order to circumvent this problem, the hyperbolic functions 
were rewritten in terms of exponentials allowing the programmed equations, in 
terms of ratios of numerator and denominator terms, to be evaluated without 
overflow problems . 

Figures 3b and 3c show the stiffness parameters needed in the theoretical 
solution for the homogeneous (i.e. solid) plate and the honeycomb plate. For 
the honeycomb plate, two different core stiffnesses were investigated. The 
stiffer one is representative of aluminum honeycomb construction that has been 
used at the GSFC. The more flexible one was chosen because it represents a 
core flexibility that is quite low and was expected to be a more critical test 
of the QUAD4's shear flexibility formulation. 

The results of the mesh study, showing the convergence of the QUAD4 solutions 
to the theoretical, are presented in tabular form in tables 1-2 and in 
graphical form in figures 4-7. Both formats show % error in displacement at 
the center of the plate as a function of mesh refinement. Results are 
included for COSMIC 88, MSC 65C and UAI 10.0 NASTRAN. In the tables results 
for COSMIC version 87 is also indicated as will be discussed below. The 
tables merely give exact numbers (along with the theoretical displacements) 
and the figures contain the same error information, but in graphic form. 

Figures 4 and 5 and table 1 are the results for the homogeneous plate. The 
difference between the results in figures 4 and 5 (and that in the two parts 
of table 1) is that figure 4 (and the top half of table 1) is for a solution 
in which shear flexibility is included and figure 5 (and the bottom half of 
table 1) is without shear flexibility. These two situations were investigated 
to test the MID3 option on the PSHELL NASTRAN bulk data deck card which allows 
the effects of shear flexibility to be ignored if MID3 is left blank. As seen 
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in figures 4 and 5 the NASTRAN results converge very rapidly with mesh 
refinement for COSMIC 88, MSC 65C and UAI 10.0. Table 1 contains the same 
information along with results for COSMIC 87, the first COSMIC version to 
contain the QUAD4 element. As seen, all versions converge to less than 0.5% 
error for a mesh size of 8 x 8 with the results without shear flexibility 
converging a little more rapidly. 

Figures 6 and 7 and table 2 are the results for the honeycomb plate. Figure 6 
(and the top half of table 2) are for the honeycomb plate with the stiffer 
core and figure 7 (and the bottom half of table 2) are for the more flexible 
core. As seen in figures 6 and 7 the NASTRAN results for COSMIC 88 and the 
two commercial NASTRAN versions converge very rapidly for the two honeycomb 
plates as they did for the homogenous plate. Table 2 contains the same 
information along with the results for COSMIC 87. As indicated, the errors in 
the first version containing the QUAD4 were extremely large for the honeycomb 
plate but, as reported above, were quite good for the homogenous plate. When 
this was discovered it was immediately reported to COSMIC. They found the 
problem in a program controlled adjustable parameter (which is used to avoid 
the infamous shear locking phenomena in earlier thick shell finite elements 
based on Reissner plate theory) and sent us a fix within two days. After 
modifying the subroutine containing the error, the results became that which 
is reported under the COSMIC 88 heading (the same fix was included by COSMIC 
in the 88 release) . 

In order to test the QUAD4's sensitivity to aspect ratio, the model with a 
12 x 12 mesh was run in which the plate side dimension in the x direction was 
varied. This causes the element aspect ratio to vary while maintaining a 
constant mesh in an attempt to remove mesh refinement errors from 
significantly affecting the results. As seen in tables 1 and 2, the QUAD4 
results with a 12 x 12 mesh (and aspect ratio of 1.0) have very little error. 
The results of the aspect ratio study are presented in figures 8-10 and 
tables 3-5. Tables 3-5 give % error in the displacement at the center of 
the plate versus aspect ratio for a model with a mesh of 12 x 12 QUAD4 
elements (over one quarter of the plate). As mentioned above, the aspect 
ratio was varied by changing the dimension of the plate along the x axis. 

Thus, the results for the aspect ratio of 10 are for a plate (and all QUAD4 
elements) that is 10 times as long in the x direction as in the y direction. 
Due to this the theoretical solution changes with aspect ratio. Figure 8 and 
table 3 are for the homogenous plate (with transverse shear flexibility) while 
figure 9 and table 4 are for the stiff core honeycomb plate and figure 10 and 
table 5 are for the more flexible core honeycomb plate. Investigation of the 
% error in the tables, as well as in figures 8-10 show that the QUAD4 has 
essentially no aspect ratio sensitivity over the range investigated. 

Based on the above results, the COSMIC QUAD4 element is seen to give very 
accurate results for the displacements in the problem investigated, both in 
comparison to the exact theory and in comparison to the two commercial 
versions of NASTRAN that we have at the GSFC. Although the results are not 
presented herein, similarity accurate results were obtained for the shear and 
moment stress resultants as well. 
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RESULTS OF TESTING USING THE MSC ELEMENT TEST LIBRARY 


As mentioned earlier, the mesh and aspect ratio studies, while a very useful 
tool in the evaluation of an element, do not test all of the important 
variables that affect accuracy in a finite element solution. The MSC element 
test library mentioned above represents a rather exhaustive series of tests 
that include many of the element related parameters which affect the accuracy 
of a finite element solution. Reference 3 gives a detailed description of the 
test problems along with theoretical answers and the results of the testing on 
several MSC elements. The reader should consult reference 3 for a complete 
description of the various problems in the test series. The portion of this 
series of element tests that relate to the QUAD4 element was run by the 
authors on the QUAD4 elements contained in COSMIC 88, UAI 9.8+ (not version 
10.0 as for the mesh and aspect ratio study) and MSC 65C. As the MSC does in 
their report, the results are presented in detail and also in a summary form 
in which the element is given a letter grade of A through F based on the 
magnitude of the error. Table 6 shows the summary results for the 15 tests in 
the series ranging from a simple patch test to modelling of beams (using the 
QUAD4 element through the depth) and various plates and shells. The meaning 
of the letter grades is given at the bottom of the table. As pointed out in 
reference 3, a failing grade for an element in one test is not a reason to 
dismiss the element. For one thing, the test scores would improve with mesh 
refinement; the mesh used in most of the problems was quite coarse. Of 
importance in this discussion is not the actual grades listed in table 6 but 
the comparison of the COSMIC grades with those from the other two programs. 

As seen in table 6, the COSMIC QUAD4 element is as good as, or better than, 
those of the commercial programs. Although not shown in table 6, the old 
QUAD2 element (included in reference 3) has a D or F grade in 9 of the 15 
problems. This is the reason for the longstanding need for an improved shell 
element and the QUAD4 element added to COSMIC NASTRAN clearly fills that need. 
Detailed results for each of the problems in the test series are contained in 
tables 7-12 and are included for completeness. 


CONCLUSIONS 


The COSMIC QUAD4 general purpose flat shell element has been shown to be an 
excellent element and significantly enhances the usefulness of COSMIC NASTRAN. 
The element has been shown to compare excellently with those available in two 
commercial versions of NASTRAN that are currently being used at the GSFC. The 
addition of an improved triangular shell element, anticipated in the near 
future, is highly desireable as a companion element to the QUAD4 in general 
analyses of complicated shell like structures. 
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List of Symbols 


a,b 

t 

P 

D, Cn.Cs 
E 

v 

G c 

AR,ARe 

w 

N X> Ny 


= plate dimensions 
= plate thickness 
= pressure load 

= plate rigidities (see Figures 3b, 3c) 

= Young's Modulus 
= Poisson's Ratio 

honeycomb core shear modulus 

aspect ratio (ratio of planar dimensions of plate or element) 
plate displacement 

number of elements in model of plate in x, y directions respectively 
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TABLE 1 

MESH STUDY 

THICK HOMOGENEOUS PLATE 
ELEMENT ASPECT RATIO 1.0 

Theoretical Displacements 
With Transverse Shear Flexibility 1 3.571x10'^ m 

(1.406xl0'3 in.) 

Without Transverse Shear Flexibility- 3.529x1 0'^ m 

(1.390x1 0‘ 3 in.) 


Mesh 

Cosmic 

87 

% Error 

Cosmic UAI 

88 Ver. 10.0 

MSC 
Ver. 65C 

With Transverse Shear Flexibility 



lxl 

12.03 

12.03 

12.03 

21.76 

2x2 

4.35 

4.34 

4.35 

2.54 

4x4 

1.67 

1.67 

1.67 

1.39 

8x8 

0.59 

0.60 

0.59 

0.53 

12x12 

0.39 

0.41 

0.39 

0.36 

Without Transverse Shear Flexibility 



lxl 

16.90 

16.83 

16.90 

26.31 

2x2 

1.12 

1.10 

1.12 

1.67 

4x4 

0.19 

0.18 

0.19 

0.50 

8x8 

0.03 

0.00 

0.03 

0.30 

12x12 

0.00 

0.03 

0.00 

0.18 
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TABLE 2 
MESH STUDY 

THICK HONEYCOMB PLATE 
ELEMENT ASPECT RATIO 1.0 

Theoretical Displacements 
G z = 1.517xl0 8 N/m2 : 2.422x1 O' 3 m 

(9.535x10*2 in.) 

G z = 1.379xl0 7 N/ m 2 : 3.102xl0' 3 m 

(1.221xKH in.) 

% Error 

Cosmic Cosmic UAI MSC 

Mesh 87 88 Ver. 10.0 Ver. 65C 


G z = 1.517xl0 8 N/m2 (22000. psi) 


lxl 

747.3 

-16.31 

-7.21 

-17.98 

2x2 

589.9 

-1.17 

4.87 

3.26 

4x4 

311.4 

-0.25 

1.46 

1.19 

8x8 

103.3 

-0.06 

0.37 

0.31 

12x12 

47.9 

-0.03 

0.16 

0.14 

■z = 1.379xl0 7 N/ m 2 (2000. psi) 
lxl -6550.4 -6.71 

10.31 

4.92 

2x2 

-5127.3 

0.26 

5.51 

4.57 

4x4 

-2689.0 

0.09 

1.42 

1.22 

8x8 

-888.5 

0.02 

0.36 

0.31 

12x12 

-412.2 

0.01 

0.16 

0.14 
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TABLE 3 

ASPECT RATIO STUDY 
THICK HOMOGENEOUS PLATE 
WITH TRANSVERSE SHEAR FLEXIBILTY 
12X12 MESH 


AR 

theoretical w, 
m (in.) 

Cosmic 

88 

% Error 
UAI 

Ver. 10.0 

MSC 
Ver. 65C 

1 

3.571xl0~ 5 

(1.406xl0- 3 ) 

-0.38 

0.39 

0.36 

2 

8.865xl0- 5 

(3.490xl0- 3 ) 

0.28 

0.26 

0.27 

5 

11.34xl0- 5 

(4.465xl0- 3 ) 

-0.83 

-0.01 

0.05 

10 

11.38x10-5 

(4.482x10-3) 

-0.04 

-0.06 

-0.02 
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TABLE 4 

ASPECT RATIO STUDY 

THICK HONEYCOMB PLATE, Gz=1.517 x 10 8 N/m 2 (22000. psi) 

12X12 MESH 


AR 

theoretical w, 
m (in.) 

Cosmic 

88 

% Error 
UAI 

Ver. 10.0 

MSC 
Ver. 65C 

1 

2.422xl0' 3 

(9.535xl0- 2 ) 

0.02 

-0.16 

-0.14 

2 

5.974x10-3 

(2.352x10-!) 

0.05 

-0.12 

-0.13 

5 

7.631x10-3 
(3. 004x10-!) 

0.24 

0.13 

0.07 

10 

7.660x10-3 

(3.016x10-!) 

0.27 

0.17 

0.14 
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TABLE 5 

ASPECT RATIO STUDY 

THICK HONEYCOMB PLATE, Gz=1.379 x 10 7 N/m 2 (2000. psi) 

12X12 MESH 


AR 

theoretical w, 
m (in.) 

Cosmic 

88 

% Error 
UAI 

Ver. 10.0 

MSC 
Ver. 65C 

1 

3.102xl0'3 

(1.221X10' 1 ) 

-0.01 

-0.16 

-0.49 

2 

7.026xl0- 3 

(2.766X10' 1 ) 

0.03 

-0.12 

0.23 

5 

8.785xlO- 3 

(3.459X10' 1 ) 

0.20 

0.01 

0.06 

10 

8.815x10-3 

(3.470X10' 1 ) 

0.24 

0.41 

0.14 
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TABLE 6 

SUMMARY OF TEST RESULTS FOR QUAD4 SHELL ELEMENTS 



Elem. . 

Loading 

IH 




Test 

In 

Plane 


KRS 

COSMIC 

88 

UAI 

9.8+ 

MSC 

65C 

1 . Patch Test 

X 


Irregular 

A 

A 

A 

2. Patch Test 


X 

Irregular 

A 

A 

A 

3. Straight Beam, Extension 

X 


All 

A 

A 

A 

4. Straight Beam, Bending 

X 


Regular 

B 

B 

B 

5. Straight Beam, Bending 

X 


Irregular 

F 

F 

F 

6. Straight Beam, Bending 


X 

Regular 

A 

A 

A 

7. Straight Beam, Bending 


X 

Irregular 

A 

A 

B 

8. Straight Beam, Twist 



All 

B 

B 

B 

9. Curved Beam 

X 


Regular 

C 

C 

C 

10. Curved Beam 


X 

Regular 

B 

B 

B 

11. Twisted Beam 

X 

X 

Regular 

A 

A 

A 

12. Rectangular Plate (N=4) 


X 

Regular 

A 

A 

B 

13. Scordelis-Lo Roof (N=4) 

X 

X 

Regular 

B 

B 

B 

14. Spherical Shell (N=8) 

X 

X 

Regular 

A 

A 

A 

15. Thick-Walled Cylinder 
(nu=.4999) 

X 


Regular 

B 

B 

F 

Number of Failed Tests (D's and F's) 


1 

1 

2 


Grading for Shell Element Test Results 


Grade 

Retirement 

A 

2% > Error 

B 

10% > Error > 2% 

C 

20% > Error > 10% 

D 

50% > Error >20% 

F 

Error > 50% 
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TABLE 7 

PATCH TEST RESULTS 


Maxium % Error in Stress 



Cosmic 

UAI 

MSC 


88 

Ver 9.8+ 

Ver. 65C 


Quad4 

Quad4 

Quad 4 

Constant-Stress Loading 

0.00 

0.00 

0.00 

Constant-Curvature Loading 

0.00 

0.00 

0.00 
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TABLE 8 

RESULTS FOR STRAIGHT CANTILEVERED BEAM 


Normalized Tip Displacement* 
in Direction of Loading 


Tip Loading 

Cosmic 

UAI 

MSC 


88 

Ver. 9.8+ 

Ver. 65C 

Direction 

Quad4 

Quad4 

Quad 4 


(a) Rectangular Elements 


Extension 

0.996 

0.996 

0.995 

In-plane Shear 

0.904 

0.904 

0.904 

Out-of-plane Shear 

0.985 

0.985 

0.986 

Twist 

0.958 

0.957 

0.941 


(b) Trapezoidal Elements 


Extension 

1.00 

0.992 

0.996 

In-plane Shear 

0.071 

0.071 

0.071 

Out-of-plane Shear 

0.980 

0.979 

0.968 

Twist 

0.937 

0.934 

0.951 


(c) Parallelogram Elements 


Extension 

0.992 

0.992 

0.996 

In-plane Shear 

0.080 

0.080 

0.080 

Out-of-plane Shear 

0.986 

0.986 

0.977 

Twist 

0.895 

0.892 

0.945 


*: Normalizing displacement values listed in Ref. 3. It is usually a theoretical value. 
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TABLE 9 

RESULTS FOR CURVED BEAM 


Normalized Tip Displacement* 
in Direction of Loading 


Tip Loading 

Cosmic 

UAI 

MSC 


88 

Ver. 9.8+ 

Ver. 65C 

Direction 

Quad4 

Quad4 

Quad 4 

In-plane Shear 

0.834 

0.833 

0.833 

Out-of-plane Shear 

0.971 

0.971 

0.951 


RESULTS FOR TWISTED BEAM 


Normalized Tip Displacement* 
in Direction of Loading 


Tip Loading 

Cosmic 

UAI 

MSC 

88 

Ver. 9.8+ 

Ver. 65C 

Direction 

Quad4 

Quad4 

Quad 4 

In-plane Shear 

0.995 

0.995 

0.993 

Out-of-plane Shear 

0.984 

0.984 

0.985 


*: Normalizing displacement values listed in Ref. 3. It is usually a theoretical value. 
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TABLE 10 

RESULTS FOR RECTANGULAR PLATE 

Normalized Lateral Deflection at Center* 



I 

Jniform Loac 


Concentrated Load 1 

Cosmic 

88 

tJAi 
V. 9.8+ 

! MSC 
V. 65C 

Cosmic 

88 

UAI 
V. 9.8+ 

MSC 

65C 

Nodes/ 

Edge 

Simple 5 

supports 

(a) Aspect 

Ratio =1.0 

2 

1.01 

1.05 

0.981 

1.05 

1.04 

1.02 

4 

1.01 

1.02 

1.00 

1.02 

1.02 

1.02 

6 

1.01 

1.01 

1.00 

1.01 

1.01 

1.01 

8 

1.00 

1.01 

1.00 

1.01 

1.01 

1.01 


(b) Aspect 

Ratio = 5.0 | 

2 

0.986 

0.983 

1.05 

0.999 

0.989 

0.811 

4 

0.988 

0.984 

0.991 

1.02 

1.01 

0.932 

6 

0.995 

0.995 

0.997 

1.03 

1.02 

0.973 

8 

0.997 

0.997 

0.998 

1.03 

1.02 

0.989 


Clamped Supports f 

(a) Aspect 

Ratio =1.0 S 

2 

1.052 

1.046 

1.008 

0.971 

0.963 

0.994 

4 

1.038 

1.034 

1.032 

1.020 

1.015 

1.010 

6 

1.024 

1.022 

1.023 

1.027 

1.018 

1.012 

8 

1.017 

1.016 

1.016 

1.013 

1.012 

1.010 


(b) Aspect Ratio = 5.0 1 

2 

1.121 

1.112 

1.314 

0.689 

0.663 

0.519 

4 

1.023 

1.019 

1.016 

0.987 

0.974 

0.863 

6 

1.013 

1.010 

1.017 

1.028 

1.019 

0.940 

8 

1.014 

1.013 

1.017 

1.034 

1.027 

0.972 


*: Normalizing displacement values listed in Ref. 3. It is usually a theoretical value. 
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TABLE 11 

RESULTS FOR THICK-WALLED CYLINDER 


Normalized Radial Displacement* 
at Inner Boundary 


Poisson’s 

Ratio 

Cosmic 

88 

Quad4 

UAI 

Ver. 9.8+ 
Quad4 

MSC 
Ver. 65C 
Quad4 

0.4900 

1.027 

1.027 

0.864 

0.4990 

1.032 

1.032 

0.359 

0.4999 

1.033 

1.033 

0.053 


*: Normalizing displacement values listed in Ref. 3. It is usually a theoretical value 
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TABLE 12 

RESULTS FOR SCORDELIS-LO ROOF 


Normalized Vertical Deflection* 
at Midpoint of Free Edge 


No. of Spaces 

Cosmic 

UAI 

MSC 

per Edge 

88 

Ver. 9.8+ 

Ver. 65C 


Quad4 

Quad4 

Quad4 

2 

1.450 

1.450 

1.376 

4 

1.070 

1.070 

1.050 

6 

1.030 

1.030 

1.018 

8 

1.019 

1.019 

1.008 

10 

1.015 

1.015 

1.004 


RESULTS FOR SPHERICAL SHELL 


Normalized Vertical Deflection* 
at Midpoint of Free Edge 


No. of Spaces 

Cosmic 

UAI 

MSC 

per Edge 

88 

Ver. 9.8+ 

Ver. 65C 


Quad4 

Quad4 

Quad4 

2 

1.020 

1.011 

0.972 

4 

1.043 

1.040 

1.024 

6 

1.023 

1.020 

1.013 

8 

1.010 

1.009 

1.005 

10 

1.004 

1.003 

1.001 

12 

1.000 

0.999 

0.998 

16 

0.998 

0.997 

0.997 


*: Normalizing displacement values listed in Ref. 3. It is usually a theoretical value. 
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Plate Size: a= 1 .0 1 6 m (40. in.)* b= 1 .0 1 6m (40. in.) 

Boundary Conditions: simply supported on all edges 
Loading: pressure load, p=6895. N/rrT 2 ( 1 .0 psi) +Z direction 
Thickness: t=0.0508 m (2.0 in.) 

*: Variable in aspect ratio studies 




AR e = a e /b e = element aspect ratio 

N x = a/2a e = number of elements in X direction in 1 /4 of plate 
N w = b/2b 0 = number of elements in Y direction in 1 /4 of plate 



Fig. 3a 


Theoretical Solution - Central Displacement 


Central Displacement 



y r 

aD L 


+ C 5 cosh(ny) + *iyC 6 sinh(py) 



sin p.x 
H 5_ 


where, 


C _ _ [ i + n 2 D(fi p — ) + 9 am tanh(am)] 

5 cosh a m L ^ V C S C n ' 2 J 

- 1 

Ce ~ 2 cosh a m 

m:c b rmc 

a m = 2 a ’ a 
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Fig. 3b 


Theoretical Solution - Homogeneous Plate Parameters 


Homogeneous Plate 
Et3 

° = 12(l- v 2) 


Cn 


5 Et 

6 v 



G= 


E 

2(1 +v) 


E = 6.89 x 1 0 1 0 N/m2 ( ] 0.0 x 1 0& lb/in2) 
v = 0.33 
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Fig. 3c 

Theoretical Solution - Honeycomb Plate Parameters 



Core Detail 


G c = 1.379 x 10^ |\j/ m 2 (2000. Ib/in2) 
or 

1.517 x 108 N/m2 (22000. Ib/in2) 
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Error in Displacement at Center of Plate, % 


Fig. 4 

Error in Displacement at Center of Plate 

Mesh Size Study 

Homogeneous Plate with Transverse Shear Flexibility 
Element Aspect Ratio 1.0 



1 2 3 4 5 6 7 8 9 10 11 12 


Mesh Size, Nx and Ny 
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Error in Displacement at Center of Plate, % 


Fig. 5 

Error in Displacement at Center of Plate 

Mesh Size Study 

Homogeneous Plate without Transverse Shear Stiffness 
Element Aspect Ratio 1 .0 
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Fig. 7 

Error in Displacement at Center of Plate 

Mesh Size Study 
Flexible Honeycomb Plate 
Element Aspect Ratio 1 .0 














Error in Displacement at Center of Plate, % 


Fig. 10 

Error in Displacement at Center of Plate 

Aspect Ratio Study 
Flexible Honeycomb Plate 
12 x 12 Mesh 


1 
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0.4 
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EIGENVALUE COMPUTATIONS WITH THE QUAD 4 CONSISTENT-MASS MATRIX 


Thomas A. Butler 
Los Alamos National Laboratory 


SUMMARY 


The NASTRAN user has the option of using either a lumped-mass 
matrix or a consistent- (coupled-) mass matrix with the QUAD 4 shell 
finite element. At the Sixteenth NASTRAN Users' Colloquium (1988), 
Melvyn Marcus and associates of the David Taylor Research Center 
summarized a study comparing the results of the QUAD4 element with 
results of other NASTRAN shell elements for a cylindrical-shell 
modal analysis. Results of this study, in which both the lumped- 
and consistent-mass matrix formulations were used, implied that the 
consistent-mass matrix yielded poor results. In an effort to 
further evaluate the consistent-mass matrix, a study was performed 
using both a cylindrical-shell geometry and a flat-plate geometry. 
Modal parameters were extracted for several modes for both 
geometries leading to some significant conclusions. First, there 
do not appear to be any fundamental errors associated with the 
consistent-mass matrix. However, its accuracy is quite different 
for the two different geometries studied. The consistent-mass 
matrix yields better results for the flat-plate geometry and the 
lumped-mass matrix seems to be the better choice for cylindrical- 
shell geometries. 


INTRODUCTION 


At the 1988 NASTRAN Users' Colloquium, results of a study 
using the QUAD4 , four-node, shell finite element for shell 
vibrations was presented (ref. 1). This study indicated that using 
the QUAD4 element with a consistent-mass matrix results in poor 
prediction of natural frequencies of a cylindrical shell. The 
errors in predicted frequencies were small for lower circum- 
ferential harmonics (n<4) and grew from approximately 10 per cent 
for the fourth circumferential harmonic to over 20 per cent for the 
sixth circumferential harmonic. The errors seemed to be relatively 
independent of the longitudinal harmonics. The authors conclude 
that the poor performance is probably caused by either a bad 
formulation of the consistent-mass matrix or, more likely, a coding 
error in the program. The QUAD 4 element is described in reference 
2 . 
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In an effort to determine whether either of the above reasons 
for the poor results is correct, a study was undertaken at Los 
Alamos to gain more insight into the problem. Earlier studies to 
evaluate the performance of the element for static problems 
indicate that the stiffness matrix formulation is correct. Also, 
results reported in reference 1 for the QUAD4 element with a 
lumped-mass matrix indicate that the problem is not with the 
stiffness matrix because the error in frequency prediction is quite 
low (less than 4 percent even for higher circumferential har- 
monics). Of course, some degradation in accuracy is expected for 
higher harmonics because the mesh density can become a limiting 
factor. 

As a first step in our evaluation, an independent check of the 
formulation and coding was performed. No problems were found with 
either the formulation or the coding. As a final check, the mass 
matrix for a single element oriented at a skew angle to the global 
coordinate system was calculated by hand, and then results of the 
code were compared. There apparently are no errors in either the 
formulation or the coding. 

A brief review of the literature on the subject of consistent- 
mass matrices does lend some insight into the problem. Clough and 
Wilson (ref. 3) state that, if the consistent-mass formulation is 
used with a displacement compatible element, resulting frequencies 
are always larger than the true frequencies. With a lumped-mass 
matrix, the frequencies may be above or below the true frequencies 
leading to the possibility that use of the lumped-mass matrix 
formulation may result in more accurate frequency predictions. The 
NASTRAN Theoretical Manual (ref. 4) describes errors associated 
with both consistent- and lumped-mass matrices for the rod 
elements. Fortunately, in this case, the errors turn out to be 
opposite in sign and of equal magnitude for lower-order terms. 
Thus, an accurate mass matrix can be generated simply by averaging 
the lumped- and consistent-mass matrices. The case does not appear 
to be the same for shell elements. Stavrinidis et al. (ref. 5) 
propose improved mass matrices for several elements, including the 
one-dimensional bar, two-dimensional membrane, and the pure bending 
beam element. Their method and other methods using so-called 
"higher-order" mass matrices depend on the predicted frequencies 
being consistently high or low. With significant effort, similar 
methods may be applicable to the current three-dimensional shell 
problem. However, as will be seen later in this paper, solutions 
with the consistent-mass matrix for the QUAD4 element can be either 
high or low, depending on the geometry of the structure. 


TEST PROBLEMS 


Two test problems were chosen for this study. The first was 
a free-free flat plate for whose natural frequencies we have 
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closed-form, analytical solutions. The second was a right, 
circular cylinder. Closed-form solutions do not exist for this 
geometry, so the finite-difference code BOSOR (ref. 6) was used 
with a fine mesh to establish the reference frequencies and mode 
shapes. BOSOR results compare favorably with approximate solutions 
presented by Blevins (ref. 7). 

The flat plate was a 10 by 10, 0.1-thick square. Its elastic 
modulus was 1.0 x 10 5 , Poisson's ratio was 0.3, and the density was 
1.0. Figure 1 shows the first three vibration modes of the plate 
with the theoretical frequencies. 

The cylindrical shell had a radius of 300 and a length of 600. 
The material thickness was 3.0. Its elastic modulus, Poisson's 
ratio, and density were 3.0 x 10 6 , 0.3, and 2.588 x 10" 4 . The 
cylinder ends were simply supported without axial constraint (rigid 
diaphragm) . Table I gives the reference frequencies calculated 
with BOSOR (ref. 6) for the cylindrical shell, along with the 
approximate solutions given by Blevins (ref. 7). 

FINITE ELEMENT MODELS 

Three different finite-element codes were used to model each 
of the two test problems. The finite-element code SPAR (ref. 8) 
was used with its E43, four-node quadrilateral element. This 
element is based on a mixed formulation first proposed by Pian 
(ref. 9). For analyzing these problems both the lumped- and 
coupled-mass matrices in the SPAR code were used. Because the E43 
element is based on an assumed-stress function, rather than an 
assumed-displacement function, its coupled-mass matrix is not 
"consistent.” That is, it is not derived from the same 
displacement functions used in deriving the stiffness matrix. Two 
types of elements were used with the ABAQUS finite-element code 
(ref. 10). The S8R5 element is an eight-node element that has only 
a consistent-mass matrix option. The S4R5 element is a four-node 
element that offers only a lumped-mass matrix. Finally, NASTRAN 
was used with the QUAD 4 element with both the lumped- and the 
consistent-mass matrix options. In addition, the problems were 
analyzed with NASTRAN using a matrix that is the average of the 
consistent- and lumped-mass matrices. 

The flat plate was modeled with three mesh densities having 
three, five, or seven nodes along each edge of the plate. ABAQUS 
was not used with the coarsest mesh because that would have 
resulted in a one-element mesh for the eight-node S8R5. The 
cylindrical shell was also modeled with three different mesh 
densities. These meshes had 5, 9, or 17 nodes on each edge. For 
the ABAQUS eight-node element, fewer total nodes were present 
because of the lack of the middle node. For this study, only one 
eighth of the shell was modeled, and symmetry conditions were used 
on all boundaries. Thus, only the even circumferential and odd 
longitudinal harmonics were determined. 
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All the NASTRAN solutions were obtained using the FEER 
eigenvalue extraction method. The mass-orthogonality test 
parameter was 0.0001 for the analyses. 


RESULTS 


Results for the flat plate are shown in figures 2 through 4. 
In these figures, the horizontal axes show the number of nodes on 
each side of the square mesh and the vertical axis is the natural 
log of the absolute value of the error in predicted frequency. The 
error is simply the ratio of the calculated frequency to the 
theoretical frequency minus 1.0. For reference, a plotted 
In (error) of -4.6 is approximately 1.0 per cent in error in 
absolute frequency determination. A plotted value of -8.0 is 
roughly equivalent to 0.03 per cent error. 

The data points labeled •'lumped,'' "consistent," and "average" 
are all for the NASTRAN QUAD4 element. Study of the results 
reveals some definite patterns. As might be expected, the 
consistent-mass matrix always outperforms the lumped-mass matrix. 
However, the rates of convergence seem to be approximately the 
same. The SPAR results that were obtained by using the coupled- 
mass matrix are consistently better than the NASTRAN results. 
However, the convergence pattern is not smooth and, for all cases, 
the SPAR E43 element with its coupled-mass matrix yields better 
answers with the intermediate, rather than the fine, mesh. This 
result is somewhat disturbing, although, in all cases, the errors 
were small. The ABAQUS S8R5 element also gives slightly better 
results than does the NASTRAN QUAD4 element. For the flat plate, 
the elements with the consistent-mass matrix formulation always 
overpredicted the frequencies and those with the lumped-mass 
matrices always underpredicted the frequencies. 

Results for the cylinder are not as clear as for the flat 
plate. Figures 5 through 7 show the frequency-convergence 
characteristics for the elements that are being considered for 
three different modes. These involve the second, fourth, and sixth 
circumferential harmonics (n=2,4,and 6) and the first longitudinal 
harmonic (m=l). The most striking observation is that, for the 
QUAD4 element, the lumped-mass matrix is now outperforming the 
consistent-mass matrix. This observation seems to confirm the 
result of Marcus (ref. 1). To illustrate the point, data from 
reference 1 have been added to the figures. Here, the definition 
of the ordinate axis has to be qualified. In reference 1, a 13 
node by 37-node mesh was used in modeling one half of a cylinder. 
This becomes a 7-node by 19-node mesh when an eighth of the 
cylinder is considered, as is the case for this study. Because, 
for the modes presented, only the first longitudinal harmonic is 
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present, we can loosely define this as a 7 by 7 mesh and plot it 
as such on our figures. 

Only for the lower, second harmonic (fig. 5) does the ABAQUS 
S8R5 element outperform the NASTRAN lumped-mass element. For all 
three modes, the QUAD 4 with the lumped-mass matrix yields the best 
results. Its deviation from the QUAD 4 with the consistent-mass 
matrix increases with higher circumferential harmonics. The SPAR, 
E43 element with a coupled-mass matrix tends to follow the QUAD4 
element closely for these modes . 

Except for a few cases, the frequencies were overpredicted 
for consistent-, lumped-, and coupled-mass matrices for the 
cylindrical-shell problem. The exceptions were the QUAD4 lumped- 
mass matrix and the E43 coupled-mass matrix for the finest mesh for 
the second and third modes considered here. 

Frequency is not the only parameter that should be considered 
for modal analyses. The other is, of course, the mode shape. One 
method of comparing mode shapes is to compare calculated 
generalized masses for the solutions using the different elements 
being considered. Another is to use a parameter frequently 
calculated when comparing calculated mode shapes with 
experimentally measured mode shapes. This parameter is called the 
mode shape correlation coefficient (MSCC) and is described in 
reference 11. It essentially provides a measure of the least- 
squares deviation of the points being compared from a straight-line 
correlation. Both these measures were used in comparing solutions 
for the n=8, m=5 mode for the cylinder being considered here. 
Results of these comparisons are summarized in table II, along with 
comparisons of the predicted frequencies. The predicted 
frequencies are normalized using the BOSOR code results as the 
baseline. The generalized masses were normalized using the 
theoretical value obtained by direct integration of the square of 
the analytically perfect mode shape multiplied by the material 
density. For the MSCC comparisons, mode shapes predicted by BOSOR 
were used as the "correct" shape. 

A study of the results summarized in table II shows again that 
the lumped-mass matrix provides better frequency predictions than 
does the consistent-mass matrix for the NASTRAN QUAD 4 shell 
element. Note that for the 9-node by 9-node mesh, the error for 
the consistent-mass matrix is over 30 per cent. A finer mesh (17 
by 17) with the consistent-mass matrix provides better frequency 
approximations, but the prediction is still not as good as for the 
lumped-mass matrix with a coarser mesh. The generalized mass is 
in considerable error for both QUAD4 cases in which the consistent- 
mass matrix is used. 

The generalized mass is a much more sensitive measure of mode 
shape error than the MSCC, as evidenced by data for the ABAQUS 
results that used the S8R5 element. Here the MSCC is quite close 
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to 1.0 for the 9-node by 9-node mesh, but the generalized mass is 
over 30 per cent in error. As the mesh is refined, the generalized 
mass improves, but it is still not as accurate as for the QUAD4 
when a lumped-mass matrix is used. Note that the performance of 
the ABAQUS S4R5 element compares favorably with the NASTRAN QUAD 4 
element . 

The SPAR E43 element, which performed nearly as well as the 
QUAD4 in predicting frequencies for all the shell modes considered 
in this study, apparently predicts both the frequency and 
generalized mass accurately if the coupled-mass matrix is used. 
However, somewhat unexpectedly, this element does not perform quite 
so well with a lumped-mass matrix. In this case, the frequency is 
predicted accurately but the mode shape has considerable error 
associated with it, as evidenced by the underprediction of the 
generalized mass. 


CONCLUSIONS 

Among the elements considered in this study, the NASTRAN QUAD 4 
element with a lumped-mass matrix seems to be the best choice when 
the geometry is a cylindrical shell. A general rule seems to be 
that, for any element considered here, consistent-mass matrices 
should be avoided for this particular geometry. On the other hand, 
for flat-plate geometries, the consistent -mass matrix outperforms 
the lumped-mass matrix. 

These conclusions imply that choices are difficult when 
modeling geometries that deviate from the simple geometries 
considered here. It is possible that an alternate method of 
deriving the mass matrix, such as the SPAR coupled-mass matrix, 
would generate a result that would be more generally applicable. 
Note that it seems to perform well for both geometries. However, 
for the present, if the geometry is predominantly cylindrical, the 
lumped-mass matrix should always be used with the NASTRAN QUAD 4 
element. 


66 



REFERENCES 


1. Marcus, M. S.; Ever stine, G. S.; and Hurwitz, M. M. : 
Experiences with the QUAD 4 Element for Shell Vibration. 
Sixteenth NASTRAN Users' Colloquium, NASA CP-2505, 

National Aeronautics and Space Administration, April 1988, 
pp. 39-43. 

2. Venkayya, V. B.; and Tischler, V. A.: QUAD4 Seminar. WRDC- 

TR— 89— 3046, Wright Research and Development Center, April 

1989. 

3. Clough, R. W.; and Wilson, E. L. : Dynamic Finite Element 

Analysis of Arbitrary Thin Shells. Computers and Structures, 
vol. 1, 1971, pp. 33-56. 

4. The NASTRAN Theoretical Manual. NASA SP-221(06), National 
Aeronautics and Space Administration, January 1981. 

5. Stavrinidis, C.; Clinckemaillie , J.; and Dubois, J: 

Concepts for Finite-Element Mass Matrix Formulations. AIAA 
Journal, vol. 27, 1989, pp 1249-1255. 

6. Bushnell , D: BOSOR 4: Program for Stress, Buckling, and 

vibration of Complex Shells of Revolution. Structural 
Mechanics Software Series - Vol. 1, N. Perrone and W. 

Pilkey, Eds., University Press of Virginia, Charlottesville, 

Virginia, 1976. 

7. Blevins, R. D. : Formulas for Natural Frequency and Mode 

Shapes. Robert E. Krieger Publishing Co., Malabar, Florida, 
1979, pp. 293-309. 

8. Whetstone, W. D. : SPAR Structural Analysis System Reference 

Manual. NASA Cr 145098-1, Engineering Information Systems, 
Inc., February 1977. 

9. Pian, T. H. H. : Derivation of Element Stiffness Matrices by 

Assumed Stress Distributions. AIAA Journal, vol 2, 1964, 
pp. 1333-1336. 

10. ABAQUS User's Manual, Version 4.7. Hibbitt, Karlsson, and 
Sorenson, Inc., 1988. 

11. Ewins, D. J.: Modal Testing; Theory and Practice. Research 

Studies Press LTD., Letchworth, Hertfordshire, England, 1986, 

pp. 222-226. 


67 



TABLE I 

PREDICTED FREQUENCIES (Hz) FOR CYLINDRICAL SHELL USING BOSOR (BLEVINS) 

FOR EVEN CIRCUMFERENTIAL HARMONICS (n) AND ODD LONGITUDINAL HARMONICS (m). 


\m 

n\ 

1 

3 

5 

2 

19.61 (21.82) 

47.97 (48.61) 

54.69 (54.83) 

4 

7.92 (8.27) 

33.27 (33.85) 

47.05 (47.30) 

6 

7.32 (7.59) 

23.64 (24.00) 

( 

39.56 (39.83) 

8 

10.76(11.68) 

20.63 (20.94) 

35.18(35.47) 


TABLE II 


COMPARISON OF FREQUENCY, GENERALIZED MASS, AND MODE SHAPE PREDICTED 
BY VARIOUS FINITE €LEMENT MODELS FOR CYLINDRICAL- SHELL MODE n=8, m=5. 


Nodes/side 

Computer code/ 
element 

Mass 

matrix 

Normalized 

frequency 

Mode shape 
correlation coef. 

Normalized 
generalized mass 

9 

SPAR/E43 

coupled 

1.032 

0.9995 

1.011 

9 

SPAR/E43 

lumped 

1.046 

0.9853 

0.810 







9 

ABAQUS/S8R5 

consistent 

0.982 

0.9998 

0.677 

17 

ABAQUS/S8R5 

consistent 

1.010 


0.971 







9 

ABAQUS/S4R5 

lumped 

0.987 


0.995 







9 

NASTRAN/QUAD4 

lumped 

1.014 

0.9991 

1.010 







9 

NASTRAN/QUAD4 

consistent 

1.336 


0.603 

17 

NASTRAN/QUAD4 

consistent 

1.066 


0.876 
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Figure 1. Mode shapes and frequencies (Hz) for flat plate. 



Figure 2. Frequency convergence for mode 1 of a flat plate 
as a function of mesh density for different finite elements. 
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Figure 3. Frequency convergence for mode 2 of a flat plate 
as a function of mesh density for different finite elements. 



Figure 4. Frequency convergence for mode 3 of a flat plate 
as a function of mesh density for different finite elements. 
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Frequency convergence for mode n=2, m=l of cylindrical 
function of mesh density for different finite elements. 
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Figure 6. Frequency convergence for mode n=4, m=l of a cylindrical 
shell as a function of mesh density for different finite elements. 
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Frequency convergence for mode n=6, n=l of a cylindrical 
function of mesh density for different finite elements. 
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NASTRAN MIGRATION TO UNIX tf \ 

by 

Gordon C. Chan 
Horace Q. Turner 
UNISYS Corporation 
Huntsville, Alabama 


ABSTRACT 


OOSMIC/NASTRAN, as it is supported and maintained by COSMIC, runs on four 
main-frame computers - CDC, VAX, IBM and UNIVAC. OOSMIC/NASTRAN on other 
computers, such as CRAY, AMDAHL, PRIME, CONVEX, etc. , is available commercially 
from a number of third party organizations. All these computers, with their own 
one-of-a-kind operating systems, mate NASTRAN machine dependent. The 30 b control 
language (JCL) , the file management, and the program execution procedure of these 
computers are vastly different, although 95 percent of NASTRAN source code was 
written in standard ANSI FORTRAN 77. 

The advantage of the UNIX operating system is that it has no machine 
boundary. UNIX is becoming widely used in many workstations, mini's, super-PC's, 
and even some main-frame conputers. NASTRAN for the UNIX operating system is 
definitely the way to go in the future, and makes NASTRAN available to a host of 
conputers, big and small. 

Since 1985, many NASTRAN improvements and enhancements were made to conform 
to the ANSI FORTRAN 77 standards. A major UNIX migration effort was incorporated 
into COSMIC NASTRAN 1990 release. As a pioneer work for the UNIX environment, a 
version of COSMIC 89 NASTRAN was officially released in October 1989 for DEC 
UITRIX VAXstation 3100 (with VMS extensions) . A COSMIC 90 NASTRAN version for DEC 
tt ttr tx DECstation 3100 (with RISC) is planned for April 1990 release. Both 
workstations are UNIX based conputers. The COSMIC 90 NASTRAN will be made 
available on a TK50 tape for the DEC UITRIX workstations. Previously in 1988, an 
88 NASTRAN version was tested successfully on a SiliconGraphics workstation. 


INTRODUCTION 


The advantage of AT&T's UNIX operating system is that it is an "open 
system", hardware independent, single and multiuser system, powerful, versatile, 
and reliable. This "open system", which may appear under different names such as 
UITRIX, XENIX, SunOS, AIX etc., is becoming the standard software today for the 
fast-growing market of workstation computers. Even IEM is going to adopt UNIX for 
its forthcoming workstations. As many more ccmputers are designed to run under the 
UNIX banner, these newcomers are getting cheaper, faster, and more powerful. The 
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result: unprecedented price competition that's making UNIX another word for cheap 
computing. The migration of NASTRAN to the UNIX "open system", is definitely the 
way to go. 


THE EARLY DEVELOPMENT 


NASTRAN is written mainly in FORTRAN language. Only about five percent of 
the source codes are machine-dependent. The early stage of migration, started in 
1984-1985, was a move towards ANSI FORTRAN 77, which is a standard FORTRAN 
compiler for all UNIX based computers. In this early stage of development, the 
NASTRAN UNIVAC version was moved from the 'FOR' compiler to the 'FIN' compiler, 
and the CDC version from FORTRAN 4 to FORTRAN 5. The VAX NASTRAN had been 
maintained as a separate version until the 1984 release. This release shared the 
machine independent source code with the other computers (IEM, CDC and UNIVAC) . 


TEST ON SiliconGraphics WORKSTATION 


A NASTRAN test program, based on OOSMIC/NASTRAN 88 VAX release, was 
converted and ran "successfully" on a SiliconGraphics workstation. Only 
occasionally this test program failed in some NASTRAN dynamic problems. Several 
UNIX job control languages (JCLs) were written to compile, link edit, and execute 
NASTRAN for this test program. These JCLs played an important part in the airrace 
of the SiliconGraphics pilot NASTRAN test. With further refinement and improvement 
(done in 1989) , the JOs, applicable to all UNIX based computers, play an 
important role in the NASTRAN migration to UNIX. The JCL to execute a NASTRAN job 
(cold start, restart or substructuring) is indeed very user friendly. 

*niis SiliconGraphics test program was also used to identify and verify 
efficiency improvements of the NASTRAN source code. The UNIX utility profiler, 
prof, was used for timing studies of the codes needing efficiency improvement. 
These studies resulted in over 30 percent speed improvement of the VAX NASTRAN 
version. The other NASTRAN versions were also benefited. 

All changes that were required to make this SiliconGraphics test program 
successful, were incorporated into the machine independent NASTRAN source codes . 


VAX NASTRAN 


The VAX version of NASTRAN is written entirely in FORTRAN language. 
Hardware-wise, VAX and many UNIX based computers are quite similar. They are 
virtual memory ccmputers with 32-bit word architecture. The file management 
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systems are quite similar. The VAX FORTRAN version is the natural choice for the 
migration of NASTRAN to the UNIX system. 

The NASTRAN GINO (General Input aNd Output file management) package of the 
VAX version has gone through extensive revision and improvement in 1987-1988. Many 
I/O processes have been shortened and streamlined. The packing and unpacking of 
matrix data were improved and speeded up. The UNIX based computers have therefore 
benefited from previous VAX improvements. (The improved VAX GINO and matrix 
packing/unpacking was also tested successfully on an IEM 3084 machine) 

The VAX 89 NASTRAN release was compiled and linked successfully on a DEC 
T TTTR TX VAXstation 3100 (with VMS extensions) , using the UNIX JCIs from the 
SiliconGraphics test program. Only one subroutine, cHJTIM that obtains the CRJ 
tiwva from the computer system, needed modification. All 119 NASTRAN demonstration 
problems, plus 20 more user problems, ran successfully. This NASTRAN UNIX version 
was officially released on a TK50 tape in October 1989. 

The DEC UITRIX DECstation 3100 (with RISC, Reduced Instruction Set Chip) 
required additional modification of the NASTRAN source codes. (See next 
paragraph. ) Occasionally this version failed in some dynamic problems, exhibiting 
the same symptom as that of the SiliconGraphics test program. There will be no 
official 89 release of this UNIX version. Presently, it is planned to have a 90 
NASTRAN release for UNIX based computers with RISC processors. 


UNIX/FORTRAN REQUIREMENTS 


The ANSI FORTRAN 77 is the standard FORTRAN compiler for all UNIX based 
computers and workstations. However, small differences may exist among ANSI 
FORTRAN 77 compilers from different manufacturers. The 1990 OOSMIC/NASTRAN 
incorporates many known specifications that are required by various ANSI FORTRAN 
77 compilers. The changes involve: 

a. External declaration of bit-shifting functions (ISHIFT and RSHIFT) , 
the logical functions (ANDF and ORF) , and complement function 
(C0MPLF) , to avoid system functions of the same names. 

b. Standardization of OPEN, READ and WRITE commands for direct-access 
files. 

c. Removal of octal and hexadecimal constants from FORTRAN executable 
source code. 

d. Elimination of jumping into an inner do-loop, which was previously 
allowed via an ASSIGN statement. 

e. Dimension of one for all open core arrays. 

f. Alignment of all open core arrays. 
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The last change (f) above is the most tricky process for UNIX based 
computers, particularly those with RISC. Throughout the NASTRAN source code, 
several hundred labelled commons are used for the open core space. (NASTRAN has no 
dimension limit. The open-ended open core is used as scratch space for internal 
computations and storage) . The other mainframe computers, particularly the non- 
virtual CDC and UNTVAC machines, require this open core space to have a unique 
name in a subroutine or group of subroutines, such that the open core space can be 
positioned strategically in the executable overlay program. To compromise among 
the virtual (UNIX based computers, VAX and IEM) and the non-virtual memory 
computers (that require program overlays), a block data routine, ZZOORE.f, was 
written to be used only for the UNIX based computers. All the open core labelled 
commons are included in this block data routine. The labelled common /XNSTRN/ must 
be the very first in the list, and /ZZZZZZ/ must be the very last. These first and 
last labelled common requirements must be true not only in the FORTRAN source 
code, but also in the compiled relocatable (or object) program. The user could use 
the UNIX command 'nm -n ZZOORE.o' to verify that /XNSTRN/ and /ZZZZZZ/ are 
positioned correctly. If they are not, something must be done to get the NASTRAN 
open core alignment correct. It is for this reason (too many labelled commons in 
one subroutine) that the DEC ULTRIX DECstation 3100 (RISC) uses two block data 
routines: ZZOORE.f for NASTRAN links 1 through 14, and ZZKORE.f for NASTRAN link 
15. It is also for this reason that a C-program, SOROBJ.c, is included in the 
UNIX NASTRAN release tape, to sort the open core labelled commons in ZZOORE.o (a 
relocatable file) , only if all other efforts fail to obtain the proper alignment. 


CONCLUSION 


The 90 OOSMIC/NASTRAN release incorporates many changes as required by the 
UNIX based computers and workstations. With a set of proven user friendly UNIX 
JCIs, it should run successfully on many UNIX based computers presently available, 
or still on the vendors' drawing boards. (FORTRAN compile and link edit are 
required.) Of course, this 90 OOSMIC/NASTRAN release will continue to operate as 
before on the IBM, VAX, CDC and UNTVAC mainframes. This is a version that bridges 
from the old world of proprietary and limited operating systems to the new UNIX 
world of "open system". 
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OBTAINING EIGENSOLIJTIONS FOR MULTIPLE FREQUENCY 

RANGES IN A SINGLE NASTRAN EXECUTION - , 

by )3 


P. R. Pamidi 
RPK Corporation 
Columbia, Maryland 


and 


W. K. Brown 
RPK Corporation 
Hayes, Virginia 


SUMMARY 

A novel and general procedure for obtaining eigenvalues and eigenvectors for multiple fre- 
quency ranges in a single NASTRAN execution is presented. The scheme is applicable to normal 
modes analyses employing the FEER and Inverse Power methods of eigenvalue extraction. The 
procedure is illustrated by examples. 


INTRODUCTION 

NASTRAN currently offers four methods for real eigenvalue extraction. They are the Tridi- 
agonal or Givens method, the Tridiagonal Reduction or FEER method, the Inverse Power method 
and the Determinant method (see Section 10 of Reference 1 for details). 

The Givens method is a transformation method that computes all of the eigenvalues in a prob- 
lem; in addition, eigenvectors corresponding to a specified range of frequencies or to a specified 
number of lowest eigenvalues can also be computed. The FEER method is also a transformation 
method that allows for the extraction of a specified number of eigensolutions. It requires the 
specification of a “shift point” or frequency around which the eigensolutions are desired. The 
Inverse Power and Determinant methods are both tracking methods that allow for the extraction of 
a specified number of eigensolutions. They both require the specification of a frequency range for 
which the eigensolutions are desired. 

When all, or almost all, of the eigenvalues in a problem are required, the Givens method is the 
generally the most efficient method to use because the total effort is not highly dependent upon the 
number of eigenvalues that are extracted. However, when the order of the problem exceeds a few 
hundred, this method may require prohibitively time-consuming out-of-core operations, thereby 
losing its efficiency. 
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If an user is not interested in obtaining all of the modes in a problem, but only in a smaller 
number of eigensolutions around a certain frequency or within a certain frequency range, the FEER, 
the Inverse Power and the Determinant methods are the obvious choices. The computations in all 
of these methods are proportional to the number of eigensolutions extracted. The FEER method is 
probably the most efficient of the three. It is quite effective in obtaining eigensolutions around the 
selected frequency. It is also very efficient even when out-of-core operations are involved. The 
results obtained by the Inverse Power and Determinant methods, on the other hand, are very 
susceptible to the number of estimated roots specified for a given frequency range (field 5 on the 
EIGR bulk data card; see Reference 2). When the specified number for the estimated roots is larger 
than the actual number of roots in that range, these methods are apt to yield many lower frequencies 
outside the specified frequency range. The Determinant method is the least efficient of all of the 
methods and will, therefore, not be considered any further for the purpose of this paper. 


CURRENT PROCEDURE FOR OBTAINING EIGENSOLUTIONS 
FOR MULTIPLE FREQUENCY RANGES 

There are practical situations in which an user may be interested in obtaining eigensolutions for 
multiple frequency ranges, with one frequency range quite distinct and apart from another frequency 
range. Complex configurations involving control systems, experimental setups and structural 
subsystems are examples of such situations. 

None of the eigenvalue extraction methods discussed earlier can accomplish the above objec- 
tive directly. So, if the user wishes to obtain eigenvalues for more than one range of frequencies in 
such cases as the above, he can accomplish it at present in one of two ways. The first way is for the 
user to make a single NASTRAN execution with a large frequency range (or a shift point in 
conjunction with a large number of desired roots) so as to encompass all of the frequencies in the 
ranges of interest. However, this will not be very cost effective if the frequency ranges of interest 
are widely separated. The alternative way is for the user to perform multiple NASTRAN executions, 
one execution for each range of frequencies, effectively utilizing the APPEND feature (see Section 
9.2.2 in Reference 1 and Section 2.3.7 in Volume 2, Reference 2). However, this latter procedure 
involves checkpoint/restart runs and is rather cumbersome for the purpose. 


PROPOSED PROCEDURE FOR OBTAINING EIGENSOLUTIONS 
FOR MULTIPLE FREQUENCY RANGES 

The above objective of obtaining eigensolutions for multiple frequency ranges can be accom- 
plished in a single NASTRAN execution by an innovative procedure that involves the use of DMAP 
ALTERs in conjunction with certain specific input data requirements. This procedure involves 
performing a normal modes analysis employing multiple subcases and using the FEER method or 
the Inverse Power method (whichever is preferred). Each subcase is setup so as to obtain 
eigensolutions in a specified frequency range. The final results of the analysis will contain the 
eigensolutions obtained in all of the specified frequency ranges. This procedure is, in essence, a 
novel application of the APPEND feature referred to above. The important difference is that, while 
the APPEND feature was originally conceived to be employed in a checkpoint/restart environment, 
the proposed procedure accomplishes this in a single NASTRAN execution. 
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The DMAP ALTER package required for the above procedure is given in Appendix A. The 
details of the input data requirements and the output obtained from the analysis are discussed below. 

Executive Control Deck 

The user should employ the DMAP ALTER package given in Appendix A either by explicitly 
including it or by referencing it via a READFILE card in the Executive Control Deck. Note that, 
for every additional subcase beyond the first subcase, the ALTER package involves a pair of OFP 
and READ modules. 

Case Control Deck 


The user should have as many subcases in the Case Control Deck as the number of distinct 
frequency ranges for which he wishes to obtain eigensolutions. Each subcase must have a separate 
METHOD request. The METHOD request in each subcase references a distinct EIGR card in the 
Bulk Data Deck that either implies (in the case of the PEER method) or defines (in the case of the 
Inverse Power method) a distinct frequency range. 


All output requests and constraint specifications must be above the subcase level. Thus, the only 
difference between one subcase and another must be the different METHOD that they request. Also, 
since the final results include eigensolutions for all of the subcases, PLOT requests should not make 
explicit references to subcase numbers. 

Bulk Data Deck 


In addition to the required modeling (geometry, constraints, etc.) data, the Bulk Data Deck 
should have as many EIGR cards as the number of subcases employed (and the corresponding 
METHOD requests) in the Case Control setup. 

When using the FEER method, the EIGR bulk data card for each METHOD request requires 
the specification of a shift point or frequency that indicates the center of a frequency range as well 
as the number of desired roots (see Reference 2). The user should specify appropriate values 
accordingly. The shift point specified has a significant effect on the actual eigensolutions extracted. 
Accordingly, depending upon the shift point specified, the actual number of roots computed may be 
more or less than the number of desired roots specified in the data. 

When using the Inverse Power method, the EIGR bulk data card for each METHOD request 
requires the specification of a frequency range, the number of desired roots as well as the number 
of estimated roots in the specified frequency range (see Reference 2). The number of estimated roots 
specified has a significant effect on the actual eigensolutions extracted. However, the user, in 
general, will not have an a priori idea of the actual number of eigenvalues that may exist in a 
particular frequency range. Accordingly, the user should use his best judgment to specify this 
number. A number for the estimated roots that is larger than the actual number of roots in that range 
will, in general, yield a number of lower frequencies that are outside the specified frequency range. 

It should also be noted here the eigensolutions resulting from any particular subcase will include 
not only the eigensolutions that are computed in that subcase, but also the eigensolutions resulting 
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from all previous subcases. Accordingly, regardless of which eigenvalue extraction method is 
employed, the number of desired roots specified on the EIGR card for a particular subcase must 
allow not only for the roots that will be computed by that subcase, but also include the roots computed 
by all of the previous subcases. 

Output from the Analysis 

As mentioned earlier, the final results from the analysis will include the eigensolutions for all 
of the subcases specified in the Case Control Deck setup. The DMAP ALTER package given in 
Appendix A also generates the eigenvalue summary table and the eigenvalues for all of the subcases 
of the analysis. If the user so desires, he can suppress any of these intermediate results by 
commenting out the OFP modules corresponding to those subcases (see Appendix A). Also, for 
every subcase, the program indicates the number of roots from all previous subcases that are included 
in that subcase. 


EXAMPLES 

Two examples were set up to illustrate the procedure discussed above. Details are given below. 
All of the runs were made on RPK’s CRAY version of NASTRAN. 

Example 1 

The standard NASTRAN Demonstration Problem No. D10-02-1A (see Reference 3) was 
selected for this example. This problem employs the Givens method for eigenvalue extraction. The 
cyclic frequencies obtained for this case are presented in Table 1. 

This problem was then modified to use the Inverse Power method and eigensolutions for two 
frequency ranges (500.0 - 1000.0 hertz and 20000.0 - 30000.0 hertz) were requested using the 
procedure described above. The input data setup is given in Appendix B. 

The cyclic frequencies obtained for this case are presented in Table 2. It can be seen that these 
frequencies are subsets of those in Table 1 . 


Example 2 

A variation of NASTRAN Demonstration Problem No. D03-08-1A (without any SUPORT 
data) (see Reference 3) was selected for this example. This problem also employs the Givens method 
for eigenvalue extraction. The cyclic frequencies obtained for this case are presented in Table 3. 

This problem was then modified to use the FEER method and eigensolutions around three shift 
points or frequencies (100.0 hertz, 1500.0 hertz and 5300.0 hertz) were requested using the 
procedure described above. The input data setup is given in Appendix C. 

The cyclic frequencies obtained for this case are presented in Table 4. It can be seen that these 
frequencies are subsets of those in Table 3. 
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The results of the above examples clearly demonstrate the validity and usefulness of the pro- 
posed method. 
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CONCLUDING REMARKS 

A novel procedure for obtaining eigenvalues and eigenvectors for multiple frequency ranges in 
a single NASTRAN execution is presented. The scheme is applicable to normal modes analyses 
employing the FEER and Inverse Power methods of eigenvalue extraction. The procedure is 
illustrated by examples. The procedure should be particularly helpful in large problems with widely 
separated frequency ranges. 
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APPENDIX A 


DMAP ALTERS for Obtaining Eigensolutions for Multiple 
Frequency Ranges in a Single NASTRAN Execution 


$ * * * * * * * * * * * ** * * * % % * * * % * * * * * * * * * * * * * * * * * * * % * % * * * * * * * * * * * * * * * * * * * * * 

$ THE FOLLOWING ALTERS ARE FOR DISPLACEMENT RIGID FORMAT 3 
$ (NORMAL MODES ANALYSIS). SIMILAR ALTERS WILL WORK FOR 
$ OTHER RIGID FORMATS THAT INVOLVE REAL EIGENVALUE 
$ EXTRACTION. 

$ 

$ NOTE THAT, FOR EVERY SUBCASE BEYOND THE FIRST SUBCASE, THE 
$ ALTERS BELOW INVOLVE A PAIR OF OFP AND READ MODULES. 

J************ *********** jIc***************.^*^^^^^^^^^^^^ 

$ 

$ INSERT AFTER THE READ MODULE IN THE RIGID FORMAT 
$ 

INSERT READ $ ON RPK-SUPPORTED VERSIONS 

$ USE ALTER 70 $ ON 1989 COSMIC-SUPPORTED VERSIONS 

$ USE OF ALTER 70 $ IS ALSO PERMITTED ON RPK-SUPPORTED VERSIONS 

$ 

$ USE THE FOLLOWING OFP STATEMENT TO REQUEST THE EIGENVALUE 
$ SUMMARY TABLE AND THE EIGENVALUES THAT ARE AUTOMATICALLY 
$ OBTAINED BY THE RIGID FORMAT FOR THE FIRST SUBCASE 
$ 

OFP OEIGS,LAMA„„//S,N,CARDNO $ 

$ 

$ COMPUTE THE EIGENSOLUTIONS FOR THE SECOND SUBCASE 
$ (THE LAST PARAMETER 2 IN THE FOLLOWING READ MODULE 
$ INDICATES THAT THE ANALYSIS IS FOR THE SECOND SUBCASE) 

$ 

READ KAA,MAA,MR,DM,EED,USET,CASECC/LAMA,PHIA,MI,OEIGS/ 
*MODES*/S,N,NEIGV/2 $ 

$ 

$ USE THE FOLLOWING OFP STATEMENT TO REQUEST THE EIGENVALUE 
$ SUMMARY TABLE AND THE EIGENVALUES FOR THE SECOND SUBCASE 
$ (THE RESULTS WILL INCLUDE THE RESULTS FOR THE FIRST SUBCASE) 

$ 

OFP OEIGS,LAMA„„//S,N,CARDNO $ 

$ 

$ COMPUTE THE EIGENSOLUTIONS FOR THE THIRD SUBCASE 
$ (THE LAST PARAMETER 3 IN THE FOLLOWING READ MODULE 
$ INDICATES THAT THE ANALYSIS IS FOR THE THIRD SUBCASE) 

$ 
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APPENDIX A 
(Continued) 


READ KAA,MAA,MR,DM,EED,USET,CASECC/LAMA,PHIA,MI,OEIGS/ 

*MODES */S ,N ,NEIG V/3 $ 

$ 

$ USE THE FOLLOWING OFP STATEMENT TO REQUEST THE EIGENVALUE 
$ SUMMARY TABLE AND THE EIGENVALUES FOR THE THIRD SUBCASE 
$ (THE RESULTS WILL INCLUDE THE RESULTS FOR THE FIRST AND 
$ SECOND SUBCASES) 

S 

OFP OEIGS,LAMA„„//S,N,CARDNO $ 


$ COMPUTE THE EIGENSOLUTIONS FOR THE LAST SUBCASE 
$ (THE LAST PARAMETER n IN THE FOLLOWING READ MODULE 
$ SHOULD BE AN INTEGER VALUE CORRESPONDING TO THE LAST 
$ SUBCASE, INDICATING THAT THE ANALYSIS IS FOR THE LAST SUBCASE) 

S 

$ THE FINAL RESULTS (WHICH INCLUDE THE RESULTS FOR ALL OF THE 
$ SUBCASES) ARE AUTOMATICALLY OUTPUT BY THE RIGID FORMAT 

S 

READ KAA,MAA,MR,DM,EED,USET,CASECC/LAMA,PHIA,MI,OEIGS/ 

*MODES */S ,N ,NEIG V/n $ 

$ 

CASE C A SECC JCA SEXX/*TR A N * $ 

EQUIV CASEXX,CASECC $ 

ENDALTER $ 
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APPENDIX B 


Input Data Setup for NASTRAN Demonstration Problem 
No. D 10-02- 1 A Modified to Use the Procedure Described 

in the Paper 


ID ... 


READFILE alters 


CEND 


SUBCASE 10 

METHOD = 100 
SUBCASE 20 

METHOD = 200 


BEGIN BULK 


EIGR, 100, INV ,500.0, 1 000.0, 5,5, ,,+EIG 1 
+EIG1,MAX 

EIGR,200,INV,20000.0,30000.0, 10, 10„,+EIG2 
+EIG2,MAX 


ENDDATA 
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APPENDIX C 


Input Data Setup for NASTRAN Demonstration Problem 
No. D03-08-1A (Without Any SUPORT Data) Modified to 
Use the Procedure Described in the Paper 


ID ... 


READFILE alters 


CEND 


SUBCASE 1 1 

METHOD =1001 
SUBCASE 21 

METHOD = 2001 
SUBCASE 31 

METHOD = 3001 


BEGIN BULK 


EIGR, 1001 ,INV, 100.0,,,5,,,+EIG 1 
+EIG1,MAX 

EIGR,2001,INV, 1500.0,,, 10, „+EIG2 
+EIG2,MAX 

EIGR,3001,INV,5 300.0,,, 12,„+EIG3 
+EIG3,MAX 


ENDDATA 
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TABLE 1 


Natural Frequencies for NASTRAN 
Demonstration No. D10-02-1A 


Mode No. 

Cyclic Frequency (Hz) 

1 

3.996 147E+01 

2 

2.364137E+02 

3 

2.504423E+02 

4 

7.01401 1E+02 

5 

7.034198E+02 

6 

1.153105E+03 

7 

1.375429E+03 

8 

1.574398E+03 

9 

1.956923E+03 

10 

2.277239E+03 

11 

2.291 26 1E+03 

12 

2.569 183E+03 

13 

2.783842E+03 

14 

2.929954E+03 

15 

3.00392 1E+03 

16 

3.41 1562E+03 

17 

4.786588E+03 

18 

6.412708E+03 

19 

8.291552E+03 

20 

1.030759E+04 

21 

1.371895E+04 

22 

1.657422E+04 

23 

2.0091 8 1E+04 

24 

2.424947E+04 

25 

2.913932E+04 

26 

3.485085E+04 

27 

4. 1 36255E+04 

28 

4.82888 1E+04 

29 

5.437856E+04 

30 

6.805413E+04 
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TABLE 2 


Natural Frequencies for NASTRAN Demonstration 
No. D 10-02-1 A Modified to Use the Procedure 
Described in the Paper 


Computed 

Mode 

No. 

Corresponding 
Mode No. 
in Table 1 

Cyclic 

Frequency 

(Hz) 

Subcase in Which 
Computed 
(See Appendix B) 

1 

2 

2.364137E+02 

First 

2 

3 

2.504423E+02 

First 

3 

4 

7.01401 1E+02 

First 

4 

5 

7.034 198E+02 

First 

5 

22 

1.657422E+04 

Second 

6 

23 

2.0091 8 1E+04 

Second 

7 

24 

2.424947E+04 

Second 

8 

25 

2.913932E+04 

Second 

9 

26 

3.485085E+04 

Second 
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TABLE 3 


Natural Frequencies for NASTRAN 
Demonstration No. D03-08-1A 
(Without Any SIJPORT Data) 


Mode No. 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 


Cyclic Frequency (Hz) 

7.341481E-04 

4. 168523E-04 

6.790643E-05 

1.5057 14E-04 

2.76597 IE-04 

6.434454E-04 

2.987344E+00 

3.372945E+00 

2.447569E+01 

2.6822 17E+01 

6.154903E+01 

7.034309E+01 

1 . 1 33579E+02 

1.1 7453 1E+02 

1. 64603 7E+02 

2.902883E+02 

2.905903E+02 

4.508246E+02 

4.515298E+02 

6.945689E+02 

6.953026E+02 

8.685333E+02 

9.752645E+02 

9.752699E+02 

1.343024E+03 

1.344180E+03 

1.466313E+03 

1.569195E+03 

1.912914E+03 

1.918371E+03 

2.02598 1E+03 

2.4465 37E+03 

2.446840E+03 

2.458828E+03 

2.742675E+03 

2.971822E+03 

3.918531E+03 

3.918534E+03 

4.451 85 1E+03 

6.261906E+03 

8.528055E+03 
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TABLE 4 


Natural Frequencies for NASTRAN Demonstration 
No. D03-08-1A (Without Any SUPORT DATA) Modified 
to Use the Procedure Described in the Paper 


Computed 

Mode 

No. 

Corresponding 
Mode No. 
in Table 3 

Cyclic 
Frequency 
(Hz) ' 

Subcase in Which 
Computed 
(See Appendix C) 

1 

9 

2.447569E+01 

First 

2 

10 

2.6822 17E+01 

First 

3 

11 

6.154903E+01 

First 

4 

12 

7.034309E+01 

First 

5 

13 

1.133579E+02 

First 

6 

14 

1.1 7453 1E+02 

First 

7 

25 

1.343024E+03 


8 

26 

1.344180E+03 


9 

27 

1.466313E+03 

Second 

10 

28 

1.569195E+03 

Second 

11 

39 

4.451 85 1E+03 

Third 

12 

40 

6.261906E+03 

Third 
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RANDOM VIBRATION ANALYSIS 

OF SPACE FLIGHT HARDWARE USING NASTRAN 

S. K. Thampi and S. N. Vidyasagar 
GE Government Services 
1 050 Bay Area Blvd. 

Houston, TX 77058 


During liftoff and ascent flight phases, the Space Transportation System (STS) 
and payloads are exposed to the random acoustic environment produced by 
engine exhaust plumes and aerodynamic disturbances. The analysis of 
payloads for randomly fluctuating loads is usually carried out using the Miles' 
relationship. This approximation technique computes an equivalent load factor 
as a function of the natural frequency of the structure, the power spectral density 
of the excitation, and the magnification factor at resonance. Due to the 
assumptions inherent in Miles’ equation, random load factors are often over- 
estimated by this approach. In such cases, the estimates can be refined using 
alternate techniques such as time domain simulations or frequency domain 
spectral analysis. This paper describes the use of NASTRAN to compute more 
realistic random load factors through spectral analysis. The procedure is 
illustrated using Spacelab Life Sciences (SLS-1) payloads and certain unique 
features of this problem are described. The solutions are compared with Miles' 
results in order to establish trends at over or under prediction. 
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INTRODUCTION 


During the past decade, the U.S. Spacelab program has made significant 
contributions to the advancement of space exploration and research. The 
Spacelab is a reusable laboratory that is carried in the cargo bay of the Space 
Shuttle. Experiments in several different disciplines such as astronomy, life 
sciences, and material science are accommodated in this modular laboratory 
for various shuttle missions. The module also contains utilities, computers, and 
work areas to support the experiments. The experiment hardware is mounted in 
instrument racks located on either side of the module, in overhead lockers, or in 
the center aisle, as shown in Figure 1. 

During liftoff and ascent flight events, the Shuttle and its payload are exposed to 
the acoustic environment produced by engine exhaust plumes and 
aerodynamic disturbances. Random vibrations are created by the response of 
the module shell to the acoustic noise inside the cargo bay. The vibrations of 
the shell are transmitted through the support structures (racks, mounting frames, 
etc.) to the payload equipment. The vibration levels that the equipment has to 
withstand depend on its own dynamic characteristics and its location inside the 
Spacelab. The equipment and its structural interfaces must be analyzed for 
these random loads in order to ensure the integrity and flight worthiness of the 
system. 

The analysis of flight hardware for random loads often relies on approximate 
formulations like the Miles' relation (ref. 1) to generate limit load factors for 
structural design. Due to the assumptions inherent in Miles' equation, the 
random vibration criteria developed through this approach tend to be over- 
conservative. In such cases, the results can be refined using alternate analysis 
techniques like time domain simulations or frequency domain spectral analysis. 
This paper describes the use of NASTRAN to perform spectral analysis to 
establish more realistic design loads. The procedure is illustrated using the 
Neck Chamber Pressure System (NCPS) assembly which will be flown on the 
Spacelab Life Sciences (SLS-1) mission. 


ANALYSIS BASED ON MILES’ EQUATION 


For a lightly damped single-degree-of-freedom (SDOF) oscillator subjected to 
random excitation through its base motion, Miles' relation is used as follows. 

N r = 3 ^ f n A Q (1) 

N r is the limit random load factor (g units), f n is the resonance frequency (Hz) of 
the SDOF system, A is the power spectral density (g 2 /Hz) of base acceleration at 
the resonance frequency, and Q is the dynamic magnification factor at 
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resonance. The quantity under the square root represents the mean square 
acceleration response and the factor 3 indicates the 3-sigma probability of 
occurrence, i.e., the probability of exceeding N r is 0.26%. 

The mean square response is the area under the response spectral density 
curve and is given by 

<u 2 > = J S(co) |H(o))| 2 dco (2) 

where H(to) is the transfer function of the system, S (to) is the base acceleration 

spectral density function, and to is the frequency in radians. The derivation of 
Miles’ relation is based on the following simplifying assumptions for evaluating 
the integral in Eqn (2). 

1 ) The actual spectral density of base excitation, S (to), is a slowly varying 
function in the vicinity of resonance. It can be conveniently approximated by a 
constant or white-noise spectral density, A = S (co n ), for computing the mean 
square response. 

2) Only the excitation energy contained within the system's bandwidth is 
transmitted. The rest is filtered away by the system and does not contribute to 
the mean square response. 

These assumptions are valid in the case of lightly damped systems with 
damping factor % « 1 . For such systems, the function |H(to)| 2 is very sharply 

peaked at to = co n , and reduces to half its peak value at a short distance, 2^to n , 
on either side of the peak. This distance, called the half power bandwidth, is 
very narrow for lightly damped systems. With the assumptions mentioned 
above, the integral in Eqn. 2 can be approximated by a rectangular area with 
base equal to the bandwidth and height equal to the product of the constant 
value of excitation spectral density and the peak value of transmittancy. This 
gives 


<u2> ~ § f n AQ (3) 

from which Miles' relation follows. 

In order to use Miles' relation for the analysis of flight hardware, the natural 
frequency of the equipment is first determined through analysis or test. The 
input random excitation spectrum for the equipment is then determined as a 
function of its location, its mounting configuration, and its mass. The input 
excitation spectrum has been established using data from previous flight or 
ground tests and is provided in the Spacelab Payload Accomodation Handbook 
(ref. 2). The spectral density at the resonance frequency of the component is 

found from this data. The dynamic magnification factor at resonance Q = 1/2 
is indicative of the system damping and is determined experimentally. For 
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components mounted on isolators, Q is determined from the manufacturers 
data on the isolator mounting. These values are substituted in Eqn (1) to obtain 
the design random load factor, N r . As the random vibration enyironment occurs 
simultaneously with other load environments during various mission phases, 
the estimated values of N r are combined with the appropriate quasi-static, 
thermal, pressure, and crew-induced load factors to generate design load cases 
for component analysis. 

Due to the assumptions inherent in Miles' relation and the idealization of the 
component as a single-degree-of-freedom resonator, the computed random 
load factors will be approximate. They tend to be overly conservative, 
especially when the natural frequency of the system is close to the peak 
frequency of the excitation spectrum. When the predicted random loads are 
unreasonably high, they lead to difficult design problems and alternate 
approaches are necessary to refine the random load estimates. 


ANALYSIS BASED ON SPECTRAL ANALYSIS 


The dynamic behaviour of large structural/mechanical systems can be 
adequately predicted only by multi-degree-of-freedom (MDOF) models. For 
linear MDOF systems, the dynamic characteristics are specified by a matrix of 
transfer functions, [H], whose elements H jk represent the ratio of steady-state 
response at point j to a sinusoidal excitation at point k. For displacement 
response 

[H(to>] = [-[M] co2 + i [B] o + [K]] -i (4) 

where [M], [B], and [K] represent the mass, damping, and stiffness matrices of 
the discrete model and o is the excitation frequency. 

The response of linear MDOF systems subjected to random excitation can be 
computed using well-established spectral analysis techniques. According to 
the theory of random vibrations, the response of a linear system with transfer 

function [H(co)], subjected to a stationary random load {P(t)}, is given by 

[S uu (co)] = [H(co)] [Spp(co)] [H*(to)F (5) 

where [S uo (to)] and [S pp (co)] are the matrices of response and excitation spectral 
density functions and * and T represent the complex conjugate and transpose 
operations, respectively. These matrices will have real auto-spectral density 
functions as their diagonal elements and complex cross-spectral density 
functions as their off-diagonal elements. By integrating the area under the 
response spectral density curve, the mean square response at any nodal point 
in the model can be obtained. The foregoing development for mean-square 
displacement response can be generalized to provide mean-square values for 
other response quantities such as velocity, acceleration, or stress. It is only 
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necessary to replace the transfer function matrix for displacement by the 
corresponding transfer function matrix for the desired response. The 
generalization also applies to the excitation which may be a point force, a 
loading condition (i.e., an ensemble of applied forces that are perfectly 
correlated), or enforced motion. 

The analysis of flight hardware subjected to random excitation can be carried 
out using the spectral analysis features of NASTRAN. A finite element model of 
the component is created which can accurately represent all its dominant 
modes in the excitation frequency range. The response calculations are carried 
out in two separate functional modules. First, the transfer function of the system 
corresponding to the desired response is computed in the Frequency Response 
module and then, the power spectral densities and other response statistics are 
computed in the Random Analysis module. The direct or modal superposition 
approaches can be used to perform the frequency response analysis. For each 
excitation source, p k , the nodal response, Uj, is determined at a series of user 

specified frequencies, coj. The ratio of output to input represents the transfer 

function element, Hjk(coj ). This is determined for each excitation source and the 
transfer function matrix, [H], is assembled from the results. 

In NASTRAN, random vibration analysis is treated as a data reduction 
procedure that is applied to the results of frequency response analysis. The 
inputs to the random analysis module are the frequency responses of desired 
output quantities due to different load sources and the auto- and cross-spectral 
densities of these random load sources. Each load source is referred to by a 
separate subcase in the case control deck and their spectral densities are 
specified as tabular functions of frequency in bulk data cards. If the sources are 
statistically uncorrelated, only the auto-spectral densities need be defined. The 
power spectral densities of response are calculated using Eqn (5) and the root- 
mean-square (rms) response is evaluated by numerically integrating the area 
under the spectral density curve. The results are printed and plotted for 
specified degrees of freedom of the model. 

As mentioned earlier, the random excitation applied to the structure could be a 
force, an enforced motion, or some other general form of excitation. In the case 
of Spacelab payloads, the random excitation is specified in terms of an 
acceleration spectrum applied at the structural support points. The "large mass" 
approach may be used to simulate this loading condition. This involves 
lumping a fictitious large mass, M a , at the degree of freedom in which the 
acceleration is to be enforced. An applied force equal to M a times the required 
acceleration is also prescribed for that degree of freedom. The inertia force is 
made so dominant through this operation that the resulting acceleration is very 
close to the required value. 
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ILLUSTRATIVE EXAMPLE 


The development of random vibration load factors for flight hardware is 
illustrated in this section using the Neck Chamber Pressure System assembly 
(NCPS). This Life Sciences Laboratory equipment item is used to study the 
effect of weightlessness on human cardiovascular control mechanisms and will 
be flown on a future Space Shuttle mission. The NCPS assembly is composed 
of an experiment enclosure which houses several components including a 
central processing unit, a motor control unit, two motor driven bellows, and a 
pressure gauge (Fig. 2). The assembly is mounted in an experiment rack using 
two support rails which are attached to the front and rear rack posts on either 
side. The front panel of the enclosure is bolted to the front rack post flanges at 
eight locations. The whole assembly weighs 48 pounds, and the installation kit 
including the slides, fittings, and fasteners weighs an additional five pounds. 

A finite element model of the NCPS assembly is constructed using mostly plate 
(CQUAD2 and CTRIA2) and bar (CBAR) elements. The model has a total of 206 
grid points and 192 structural elements.The masses of internal components are 
lumped at the respective centers of gravity, and stiff bar elements are used to 
connect them to the attachment points. The fasteners are modelled using rigid 
elements. Eigenanalysis was performed on the model with free boundary 
conditions to verify that the model has six rigid body modes. The analysis is 
repeated, with the rack-to-component interface points appropriately 
constrained, in order to determine the flexible modes of the component. The 
first twenty frequencies of the constrained model are shown in Table 1 . An 
inspection of the modal deformation plots and mass participation factors shows 
that the first system mode in X direction is 80 Hz. The power spectral density of 
the input excitation corresponding to this frequency is found to be 0.02 g 2 /Hz 
from Figure 3. For a conservative estimate of Q = 10 for the dynamic 
magnification factor, Miles' relation (Eqn 1) yields 

N rx = 15.04 g units (6) 

Random load factors, N^ and N rz , are computed in a similar manner. 

Random load factors can also be determined through spectral analysis. The 
computation of N rx is described here for comparison with the Miles' approach. 
The transfer function of the system is first determined by applying a unit 
sinusoidal load in the X direction at each point where the NCPS interfaces with 
the rack structure. The load is applied through the DLOAD and RLOAD cards 
which in turn refers to DELAY, DPHASE, DAREA, and TABLED cards. For a 
constant phase, unit sinusoidal input, the DELAY and DPHASE cards may be 
omitted and a unit value specified for all frequencies on the TABLED card. The 
input acceleration spectrum (Fig. 3) is specified through RANDPS and TABRND 
cards. This random acceleration is enforced at all the interface points in the X 
direction using the large mass approach. A fictitious mass, M a , about 1000 
times larger than the existing grid point mass, is lumped at the interface X 
degree of freedom using a CMASS2 card and an equal value is specified in the 
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corresponding DAREA load card. This produces the desired acceleration which 
can be verified by plotting the acceleration spectrum at the input points and 
comparing it with Figure 3. If the agreement is not adequate, the value of M a is 
increased until a good match is obtained. 

The direct and modal solution techniques are used to carry out the response 
calculations. The relative efficiency of the two approaches depends on several 
factors including the number of modes retained in the analysis and the number 
of response frequencies. The selection of these parameters always represents 
a compromise between accuracy and efficiency. The frequencies chosen for 
response computations should have good resolution in the vicinity of system 
resonances in order to obtain reliable estimates of rms response. A total of 150 
frequencies in the 0 to 500 Hz range are used for this analysis. The 
specification of damping properties in the direct and modal formulations are 
somewhat different. In the direct formulation, structural damping proportional to 
the stiffness matrix terms is specified both on the material data cards and as an 
overall uniform damping factor on the PARAM G data card. In the modal 
formulation, the damping factor is specified as a tabular function of frequency 
through the SDAMPING and TABDMP1 cards to represent the variation of 
structural damping for different modes. A damping factor of 0.1 is used for this 
analysis in order to be consistent with Q = 10 used in Miles’ relation. 

The set of output points at which the response power spectrum and rms values 
are to be recovered must be chosen judiciously. The selection can be based on 
the same criteria used for choosing an ASET (analysis set of dynamic degrees 
of freedom); namely, the points should be uniformly dispersed throughout the 
structure and should include all large mass items. A set of 12 output points 
were chosen for the NCPS and the rms acceleration response at these points is 
computed (Table 2). When these values are averaged and the 3-sigma 
probability of occurrence criteria is applied, one obtains 

Nrx = 10.67 g units (7) 

This is almost 30% less than the value predicted using Miles' relation. 

The rms response calculated using the direct and modal solution techniques is 
summarized in Table 2. For the same number of response frequencies, the 
modal solution required 350 seconds of CPU time with 20 modes included, 450 
seconds with 40 modes, 600 seconds with 60 modes, whereas the direct 
solution required 13160 seconds. The response spectrum at selected output 
points on the model is shown in figures 4, 5, and 6. The solution obtained with 
20 modes is clearly inadequate while the plots obtained using 60 modes and 
the direct approach are virtually indistinguishable. The spectra 
characteristically peak at 80 Hz which corresponds to the first X mode of 
vibration. The effects of higher modes can also be seen in the spectral plots. 
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CONCLUDING REMARKS 


An alternate method of estimating random vibration load factors for design and 
analysis of Spacelab payloads is presented. This method, based on spectral 
analysis, yields more refined random load estimates at the expense of being 
computationally more intensive than the Miles' approach. The computational 
effort can be reduced by using the modal formulation rather than the direct 
formulation for analysis. Significant reductions can be obtained in the random 
load estimates using this method. While the actual reduction depends upon the 
payload configuration being analyzed, reductions of 20 to 30% are typical. This 
method could be used to resolve difficult design problems owing to 
unreasonably high random load predictions by the Miles’ relation. 
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TABLE I. 


EIGENVALUE ANALYSIS SUMMARY 


MODE 

NO. 


EIGENVALUE 



RADIAN 

FREQUENCY 


2 . 367843E+02 
2 . 937815E+02 
3 . 146429E+02 
3 . 505287E+02 
3 . 851169E+02 
4 . 355316E+02 
4 . 797886E+02 
5. 035146E+02 
6. 316066E+02 
6 . 5042 14E+02 
6. 658782E+02 
6 . 846459E+02 
7 . 141465E+02 
7 . 699639E+02 
8 . 793395E+02 
8 . 928888E+02 
9. 324113E+02 
9 . 945816E+02 
1. 050721E+03 
1. 109605E+03 
1 . 221881E+03 


CYCLIC 

FREQUENCY 


3 . 7 68539E+01 
4 . 675677E+01 
5. 007697E+01 
5. 578837E+01 
6. 129326E+01 
6. 931701E+01 
7 . 636072E+01 
8 . 013683E+01 
1. 005233E+02 
1. 035178E+02 
1. 059778E+02 
1. 089170E+02 
1. 136599E+02 
1 . 225435E+02 
1. 399512E+02 
1. 421077E+02 
1 . 483979E+02 
1. 582926E+02 
1. 672274E+02 

1 . 7 65992E+02 
1 . 944685E+02 


99 



TABLE II RMS ACCELERATION RESPONSE SUMMARY 


GRID 

POINT 

LOCATION 

MODAL SOLUTION 

DIRECT 


20 MODES 

(g) 

40 MODES 

(g) 

60 MODES 

(g) 

SOLUTION 

(g) 

5059 

Rear Left 
Side Panel 

2 . 07 

3.50 

3.52 

3 . 52 

5080 

Rear Right 
Side Panel 

2.06 

3.97 

3 . 99 

3 . 98 

5056 

Front Left 
Side Panel 

2.47 

2.47 

2 . 49 

2 . 50 

5077 

Front Right 
Side Panel 

2 . 47 

2.47 

2.49 

2 . 50 

5076 

Top Left 
Rear Panel 

2 . 08 

4 . 53 

4.58 

4 . 58 

5073 

Top Left 
Front Panel 

2 .24 

2.48 

2.87 

2 . 88 

5094 

Top Right 
Front Panel 

2.24 

2.48 

2.87 

2 . 88 

5097 

Top Right 
Rear Panel 

2.08 

4.89 

4.91 

4.91 

5193 

M/C Unit 

2 . 07 

4.33 

4 .46 

4 . 46 

5204 

CPU Unit 

2 .81 

3.89 

3.95 

3 . 95 

5205 

M/C Mount 

2.12 

4.03 

4.04 

4.03 

5042 

Front Panel 

2.47 

2.47 

2.49 

2 . 50 

Average 


2.26 

3.46 

3 . 55 

3 . 56 
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TYPICAL SPACELAB CONFIGURATION 
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FIG. 5 . ACCELERATION RESPONSE SPECTRUM AT GRID POINT 5205 
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OBTAINING AN EQUIVALENT BEAM 

THOMAS G. BUTLER 
BUTLER ANALYSES 


In modeling a complex structure I was faced with a 
component that would have logical appeal if it were modeled as a 
beam. It was a mast of a robot controlled gantry crane. The 
structure up to this point already had a large number of degrees 
of freedom, so the idea of conserving grid points by modeling the 
mast as a beam was attractive. I decided to make a separate 
problem of the mast and model it in three dimensions with plates 
then extract the equivalent beam properties by setting up the 
loading to simulate beam like deformations and constraints. The 
results could then be used to represent the mast as a beam in the 
full model. This seemed to be a straight forward approach, but 
it was sufficiently challenging that it merited publishing a 
paper on this topic. 

The endeavor is to obtain the area A, the area moments of 
inertia II and 12, and torsional area moment of inertia J of a 
prismatic beam that would be an equivalent of the crane mast over 
its full length. The detailed model involved about 4500 uncon- 
strained degrees of freedom. The mast structure was essentially 
a hollow steel tube of square section with a cylindrical indenta- 
tion along its length on one surface only. Complications that 
made it difficult to estimate equivalent properties analytically 
were the placement of two types of interior partial shear stiff- 
eners at regular intervals along its length. These two different 
types of shear stiffeners alternated on opposite sides from each 
other most of the length. This posed no difficulty to model 
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STAINING '^EQUIVALENT BEAM 


elastically in a three dimensional model. The interesting phase 
is the loading of the 3-D model in order to simulate beam action. 

To put the problem in perspective, review for the moraen r , 
the definition of beam stiffness. 

DEFINITION: Beam stiffness is the array of forces pro- 
duced at the six degrees of freedom on both ends when a 
single degree of freedom at one end is deformed a unit 
amount while enforcing all other eleven degrees of 
freedom at both ends to be zero. 

But the Bernoulli Euler formulation of the beam as used in finite 
element analysis programs does not faithfully follow this 
prescription of stiffness to the letter.' For example, when one 
end is displaced a unit transversely, action is assumed to occur 
in - plane only. Diagrammatically the boundary conditions of the 
centroid of the B.E. formulation are indicated in the sketch. 



Note that the length remains invariant, because its transversely 
deformed end is not constrained in the axial direction. In 
effect, with this B.E. approach, the end position contracts when 
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OBTAINING AN EQUIVALENT BEAM 


bending deformation occurs. This is shown in exaggerated fashion 
in the sketch. 


Def 



If the length is not allowed to deform, Poisson deformation does 
not occur and therefore needs no constraining force to inhibit 
Poisson deformation. But if the true definition of beam stiff- 
ness were adhered to in the finite element beam, the axial posi- 
tions of the ends would be held to zero displacement and the beam 
would lengthen as transverse deformation occurs. Such axial 
stretching would result in Poisson contraction in both transverse 
directions. But if transverse translational deformations were 
held to zero, as the definition of stiffness demands, such con- 
straints would exert forces to prevent Poisson contraction. For 
instance, the transverse forces at the end of a solid beam of 
square section with a full set of constraints applied would 
appear as sketched. 

■in— Undz't' 

P^/ : 56Tr/ 


tttttt 



Dim 
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The dilemma now is to try to define what kind of equivalence 
should be sought. If A, II, 12, and J were obtained with true 
stiffness constraints, would it be proper to operate as an 
equivalent beam according to those entries on the property card, 
so as to exclude bending/axial coupling even though such action 
was present during the sample run? Or would it be more proper to 
use only B.E. conditions to get the properties that will used as 
a B.E. beam? If the latter were chosen, the question arises as 
to how faithfully we would be representing equivalence to the 
true structure. Having some doubts as to how to proceed, I 
modeled the constraints in two different ways; with full end 
constraints and with B.E. end constraints and compared the 
results. The sketch shows the constraints imposed for the two 
models. One of the things to consider in the B.E. simulation is 
that the theory requires planes to remain plane in bending. 
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The next question is: After constraint forces are 
measured, will it be acceptable to derive sectional properties by- 
substituting into the formulation based strictly on B.E.? Thar 
is to say, should the stiffness forces obtained on the left be 
equated to the B.E. formulas on the right? Just enough of the 
matrix on the right is shown to illustrate the problem. 



3 

GJ/L 


Not having any reference to use for the fully coupled 
beam I chose to use B.E. formulation to evaluate sectional 
properties for both types of modeling. 

The next question is: After accepting B.E. formulation, 
what basis should be used to reconcile differences in results of 
the methods? The reconciliation method is to use an estimation 
of the computed value of the section without the shear panels 
present as per the dimensions in the sketch. 
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COMPARISON OF PROPERTIES DERIVED FROM MODELS OF DIFFERENT 
CONSTRAINTS VS MANUAL CALCULATIONS 


SOURCE | 

A | 

11 1 

12 | 

J 

FULL CONSTRAINTS | 

36.72 | 
| 

2,444.47 | 
1 

2,605.14 | 
| 

INVALID 

1 

B.E. CONSTRAINTS | 

1 

36.66 | 

1 

728.63 | 
| 

i 

853.62 | 
| 

600.19 

1 

MANUAL | 

1 

35.96 | 

1 

2,541.82 | 

1 

3,517.62 | 

4,710.60 


This exercise had some unexpected results. The whole 
purpose of the exercise was to get an equivalent beam by using a 
full 3-D model instead of making an analytical estimate because 
of the uncertainty in being able to represent the effect of the 
partial shear panels correctly. One expects that the effect of 
the shear panels is to stiffen the steel tube, but the 3-D re- 
sults showed less stiffness than the manual check which neglected 
the panels. Why? 

In going back to examine the axial displacements in the 
3-D model using the B.E. constraints, it indicated that the end 
faces tilted instead of remaining perpendicular to the undeformed 
centroidal axis as the B.E. theory requires. The total burden of 
meeting the requirement of zero slope at the displaced end was 
put on the QUAD4 elements which formed the side panels of the 
steel tube. That is; the open ended tube had two surfaces that 
could carry such bending and two surfaces unable to carry 
in-plane shear about their normals. Even those that picked up 
such bending couldn't transmit this moment to the QUAD'S on the 
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perpendicular surface, 30 an inadequate moment developed to 
produce a net slope of zero at the centroid. By letting the end 
points displace axially they were in no position to create 
couples to satisfy the moments for zero slope. That is why the 
model with the B.E. constraints produced inadequately stiff 
sectional properties in bending and torsion. 

Going back to the model with fully constrained ends, the 
explanation as to why this model was also inadequate for simulat- 
ing an equivalent beam was this. Even though it did develop 
couples which formed the resisting moment for zero end slope by 
holding the axial displacements to zero; it still felt the defi- 
ciency of moments about the normals of side panels. In effect 
membrane action on corner displacements alone was not sufficient 
to represent the true structure without the help of the existing 
— but unrepresented — in plane shear from moments about the 
normals of the panels. 

In the case of torsion the fully _onstrained model was 
invalid because it developed local equilibrium at the end under- 
going unit rotation. The unit rotation about the axial direction 
for every end grid point was inhibited by the translational 
d.o.f.'s being held to zero. The deformation became a scalloped 
pattern instead of a uniformly rotated face. Representation of 
torsion with the B.E. model was also inadequate because it re- 
quired, but didn't get, the assistance of the panels on all four 
sides to carry the rotation about their normals. 

Does this mean that if no attempt were made to model the 
mast as an equivalent beam, but a full 3-D model were used, that 
the 3-D model would be invalid? Not at all. What it shows is 
that the 3-D model is ineffective in trying to conform to the 
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requirements of an equivalent beam representation. If a full 3-D 
plate model were used in the complete representation of the crane 
structure, good results would be obtained. 

Since the attempt is to economize on the size of the 
model, a better way to achieve the same results is to use sub- 
structuring and condense the mast to equivalent end boundary and 
intermediate mass points. 

The spirit in which this paper is presented is to publish 
failures as well as successes to help analysts avoid retracing 
the ground that has already been plowed. 
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ABSTRACT 

A nonlinear elastic force— displacement relationship is used to calculate the 
transient impact force and local deformation at the point of contact between impactor and 
target. The nonlinear analysis and transfer function capabilities of NASTRAN are used to 
define a finite element model that behaves globally linearly elastic, and locally nonlinear 
elastic to model the local contact behavior. 

Results are presented for two different structures: a uniform cylindrical rod 
impacted longitudinally; and an orthotropic plate impacted transversely. Calculated 
impact force and transient structural response of the targets are shown to compare well 
with results measured in experimental tests. 


INTRODUCTION 

Aerospace structures are subjected to impact loading from a variety of sources, 
including dropped tools, runway debris, and munitions. In some advanced materials, even 
low velocity impact can cause significant structural damage. Therefore, development of 
accurate means of calculating structural response due to impact loading can be of critical 
importance. In this paper, a computational technique is developed to predict the dynamic 
response of a structure to low velocity elastic impact. 

Structural damage due to impact invariably initiates in the immediate vicinity 
of the impact. Therefore, it is important that the local stress field in the region of contact 
be calculated accurately. Hertz [1] derived an elasticity— based force— displacement 
relationship that describes contact between two elastic bodies. The Hertzian contact law is 
given by: 


F = Ko" 


( 1 ) 


where 


F = elastic contact force 
K = contact stiffness 
n = exponent 
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and 


a = relative displacement (indentation) between impactor and target 
= Ui — Ut (i = impactor, t = target) 

The exponent n was shown in reference [1] to have the value of 3 / 2 . In dynamic 
applications such as this; F, u, and a are all time— varying. 

During low velocity impact, where impact damage is confined to the area 
immediately around the point of contact, areas of the structure remote from the impact 
may still deform in a linear elastic manner. An efficient finite element model, therefore, 
would combine a linear elastic model of the global structure with a non— linearly elastic 
behavior at the point of contact. The nonlinear force— displacement relationship in 
equation (1) is incorporated into a linear elastic finite element model (MSC/NASTRAN 
transient solution 27, COSMIC/NASTRAN transient solution 9) by using a NASTRAN 
transfer function definition and nonlinear analysis capability. In the following section, the 
Hertz contact law is discussed in addition to a method of incorporating it into NASTRAN. 
Impact loading of two different structures is then analyzed. The first problem is a 
one— dimensional rod of uniform cross section impacted longitudinally. The second is an 
orthotropic plate under transverse impact. 


CONTACT LAW 

In reference [2] Hertz derived the force— displacement relationship for two 
spherical isotropic elastic bodies of radius ri, and r 2 in contact: 

F = K a 3 / 2 (2) 


where 


is the contact stiffness and 


K - “4” 

A 

ri r 2 

k, k 2 

r 1 + r 2 
E- 

ki + k 2 

k. = - 

t 

i= 1,2 

t 

1 — iP'- 

t 

Young’s modulus and 

poisson’s 


( 3 ) 


( 4 ) 


subscripts 1 and 2 refer to each of the spheres. When a spherical impactor contacts a flat 
target, (3) simplifies to 


K 



kjkt 
ki + k t 


( 5 ) 


where i and t represent the impactor and target respectively and the k t and ki are given by: 
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( 6 ) 


kt = 


Et 

r 

1 -u t 


ki 



( 7 ) 


In equation (2), a is the local indentation at the contact point, shown 
schematically in figure 1. We have: 


a = ui-u t (8) 

where a is the relative local displacement between impactor and target at the point of 
contact. 


NASTRAN Implementation 

The non-linear local behavior was incorporated into the NASTRAN finite 
element model as follows: 

The impactor is modeled as a lumped mass just touching the target at t=0 and 
with an initial velocity towards the target. The difference between the displacement of this 
lumped mass and the displacement of the target is the indentation, o. The modeling of the 
contact between impactor and target is performed by utilizing the transfer function card, 
TF, and the nonlinear force card, NOLIN3. The TF card acts as a dynamic multipoint 
constraint, relating the displacement, velocity and acceleration of several independent 
degrees of freedom to a dependant degree of freedom. In the work discussed here, only 
displacement relationships were used. On the TF card coefficients of the following 
equation are specified [3]. 


(B„ + B,p + B,p> dep + (A' 0 + A J lP + A , ,py i „ d = 0 (9) 

j = 1 


where 


B 0 , Bj, B 2 = the coefficients for the dependant degree of freedom 

Ag, A J j, Aj = the coefficients for the independent degrees of freedom 
u dep = displacement of the dependant degree of freedom 

u ind = displacements of the independent degrees of freedom 

n = the number of independent degrees of freedom 

3 ^2 

p = the differential operator — , and p 2 = — ^ 

For this analysis, the equation would appear: 
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(l.O)ilextra point + ^( — 1-O)uiinpactor + (l.O)utargetj — 0 


( 10 ) 


that is 


n = 2 


Bp Bj, Ap A 2 
B n 


0.0 (j = 1, n) 
1.0 


Aq = -1.0 
A 0 2 = 1.0 


The resulting equation defines the indentation at every time step and assigns the value to 
an EPOINT. The EPOINT, or extra point, is used as a nonstructural variable that 
provides a place to store the value of the indentation. The EPOINT is provided as input 
to the NOLIN3 card. 


The NOLIN3 card is the means of applying the time— dependent nonlinear load 
based on the indentation. The NOLIN3 card has the form: 


where 


p(t) 


SM 4 )) -4 . > o 

0 , x(t) < 0 


( 11 ) 


P(t) = is the resulting nonlinear force 
S = is a scale factor 

x(t) = is the displacement or velocity of a degree of freedom 
A = is an amplification factor 


In modeling of the impact, we define x(t) to be the displacement of the EPOINT, S to be 
the Hertzian stiffness, and A to be 3 / 2 , as given in equation (2). Recall that the 
displacement of the EPOINT is really the indentation as obtained from the TF card. The 
resulting function then has the form: 


p(t) = 


K(a(t))’<\ o(t) > 0 
0 , a(t) < 0 


( 12 ) 
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Note that when a is less than or equal to zero (ie. the target and the impactor are out of 
contact) then the force is also zero. Two NOLIN3 cards are used, one to apply the impact 
force to the target and the other to apply the same force to the impactor in the opposite 
direction of its initial velocity. This methodology allows the impactor to slow with 
increasing impact force and eventually to unload the target as the impactor begins to travel 
in the opposite direction, away from the target. 


RESULTS 

One Dimensional Rod 

The first problem analyzed is the longitudinal impact of a steel ball on a long 
aluminum rod of constant cross section. Geometry and material properties of the impactor 
and target are given in figure 1. The problem was modeled using 144 1— D rod elements 
with each grid point having a single longitudinal degree of freedom. Two more degrees of 
freedom were used to model the impact dynamics, resulting in a total of 147 degrees of 
freedom. A single lumped mass with an initial velocity was used to represent the impactor. 
The Hertzian force— displacement relationship in equation (1) was prescribed using the 
NASTRAN NOLIN3 card, as shown in the example input deck in the appendix. 

The calculated impact force history compares well with experimentally 
determined values [4], as shown in figure 2. The calculated strain response at the midpoint 
of the target bar is compared with measured values in figure 3. The sign reversal of the 
second pulse is caused by the reflected tensile stress wave generated by the incident 
compressive wave reaching the free end of the bar [5]. 

Some insight into the timing and the location of the impact— induced structural 
failure can be gained by tracking the distribution of energy in the impactor and the target, 
as shown in figure 4. The energy balance can be expressed as: 

U t = KEi + SEi + KE t + SE t (13) 

where 

U t = total energy in system 

KEi = impactor kinetic energy = l / 2 mvj (14) 

SEi = impactor strain energy = J F(a) da = 2 /s Ka (15) 

KE t = target kinetic energy (16) 

-j (n = number of elements) 

(17) 

(n = number of elements) 


n - 1 




m J 


Vi + Vj»: 

2 ^ 


SE t = strain energy of target 


n - 1 




kj ( <5j — <5j ) 
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The total energy in the system, Ut, is divided between the kinetic energy and 
the strain energy of the target and impactor in a time— varying manner. Because damping 
effects are not considered, the total system energy is constant and equal to the initial 
kinetic energy of the impactor. The strain energy of the impactor is non— zero only during 
the contact interval (0 < tC/L < 0.4, where t = time, L = the length of the bar, and C = 
the wave speed in the bar) and peaks when the contact force is greatest, approximately 
halfway through the contact interval. The kinetic energy of the impactor decreases rapidly 
as the impactor slows during contact with the target. Eventually, at tC/L = 0.25, the 
impactor velocity (and therefore its kinetic energy) decreases to zero and the elastic 
rebound begins. The kinetic energy of the impactor never returns to its initial level 
because approximately 80% of the energy has been transferred to the target in the form of 
strain energy and kinetic energy. The strain and kinetic energies in the target both 
increase rapidly during the contact with the impactor and remain constant after contact 
has ended (tC/L > 0.4). Both strain and kinetic energies maintain equal and constant 
values until the compressive stress wave generated by the impact reaches the far end of the 
free— free bar (tC/L = 1.0). A tensile stress wave is generated when the compressive pulse 
reflects from the stress free boundary [5]. The superposition of the incident and reflected 
pulses momentarily leaves the bar stress— free which causes the strain energy to decrease to 
zero. The kinetic energy simultaneously increases, maintaining a conservation of total 
energy. The reflection process is repeated at tC/L = 2.0, when the reflected pulse returns 
to the other end of the bar. Similar energy dissipation diagrams may prove useful in 
analyzing dynamic failure of more complex structures. 


Composite Plate 

The low velocity transverse impact of a composite plate made from Scotchply 
1003 prepreg [6] is now analyzed. The problem is depicted schematically in figure 5, and is 
described in detail in references [7,8]. A modified Hertzian contact stiffness has been 
proposed [9] for application to composite materials. Specifically, equation (6) is replaced 
by 


kt = E 33t (18) 

where E 3 3 t is the transverse modulus of the plate. Plate membrane and bending stiffness 
material properties were calculated using the COBSTRAN (Composite Blade Structural 
Analyzer) computer code [10] which calculates elastic moduli of composite materials from 
known constituent properties and laminate ply orientations. 

A uniform square mesh of QUAD4 elements was used to model the 15.24 cm * 
15.24 cm (6 in * 6 in) target plate. A mesh convergence study was performed to establish 
the degree of mesh refinement necessary to arrive at a numerically converged solution. 
Three different meshes were considered, 25 * 25, 49 * 49, and 61 * 61 elements. Of these, 
the latter two produced essentially the same strain response for a given impact velocity and 
were therefore considered to be converged solutions. The results presented here were 
therefore calculated using the 49 * 49 element model. Five degrees of freedom (u x , u y , u z , 
0 X and 0 y ) were used at each nodal point, giving the model a total of 11510 degrees of 
freedom. The problem was solved on a Cray XMP in 52 CPU minutes. 
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The impactor used in the tests [7.8] was a uniform 2.54 cm (1 in) long, 
blunt— ended steel rod of radius 0.047625 cm ( 3 /ie in)- In the analysis a contact radius of 
0.047625 cm ( 3 /i 6 in) was assumed in the Hertzian contact stiffness calculations. The 
calculated impact force history is shown in figure 6. Although no direct measurement of 
the impact force was obtained experimentally, the contact time was measured [8] and found 
to be 204 microseconds. This is in good agreement with the calculated result. Figure 6 
also shows that a secondary impact occurs during the latter half of the contact interval (t 
= 175 fjs ec), probably due to the vibration of the target plate during contact with the 

projectile. 


The resulting displacement response of the plate is shown in figure 7, where it 
has been assumed that no damage occurs in the target during contact with the impactor. 
This assumption is valid based on the available test data. Ultrasonic C— scans of the 
specimens after impact indicate that this level of impactor kinetic energy (10 Joules) is 
very near the threshold energy level required to cause damage [8] in specimens of this 
layup. As a result, very little damage occurs at this impactor velocity. The anisotropic 
bending stiffness of the target (figure 5) is evident from the elliptical displacement 
contours, as the flexural disturbance travels faster in the stiffer direction (figure 7). 

The strain response at gage A is compared to the calculated response in figure 8. 
The two curves are similar in amplitude and duration but the calculated strain appears to 
lag the measured values by approximately 25 microseconds. This may be due to the 
difficulty in establishing experimentally the precise time at which contact occurs based on 
strain gage readings taken at some distance from the point of contact. The comparison 
shown in figure 9 for gage B likewise shows a time shift of approximately 25 microseconds 
between the measured and the calculated response. The amplitude and duration of the 
calculated strain response correlate quite well with the measured signal. 


SUMMARY 

A simple means of modeling low velocity, non— damaging impact using 
NASTRAN was demonstrated. A nonlinear elastic contact model was included in the finite 
element analysis using NASTRAN transfer function definitions and nonlinear analysis 
capabilities. The same contact law was used to define the force— indentation relationship 
for two different impactor/target combinations. Results in both cases showed that the 
impact force and resulting transient structural response of the target compared well with 
experimentally measured values. 
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Figure 2: Impact Force for Bar Problem 
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Figure 4: Energy Distribution for Longitudinal Bar Impact Problem 
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Figure 5: Composite Plate Impact 
Specimen Configuration 
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ANALEX - STRUCTRAL MECHANICS BRANCH 


ID TRANS. LOAD 
APP DISP 
TIME 60 
SOL 9 
CEND 

TITLE - COSMIC: TRANSIENT RESPONSE ANALYSIS: HERTZIAN IMPACT FF 


SUBTITLE- 36" AL. ROD 5/8 STEEL BALL V0-62 . 1 IN/S 

LABEL - ROO ■ • < IMPACT 

$ NONLINEAR LOAD 

NONLINEAR - 5 

$ INITIAL CONDITIONS SET 

IC - 1 

TFL-111 

SPC - 4 

TSTEP - 7 

$ OUTPUT STUFF 


SET 30 - 1.72.73.999.1001 
NLLOAD - 30 
STRESS (PR I NT) - 30 
DISP(PRINT) - 30 
BEGIN BULK 

$ * 

$ EXTRA POINT - INDENTATION 

EPOINT , 1001 

GRID. 999. .-0.3125.0.0.0.0 
GRID. 1 . .0.0, 0.0. 0.0 
-(144). .(1).-. •(0.25).— 

CROD . 1 , 1 . 1 . 2 
—(143) .— (1).— .-(1)»*(1) 

$ LUMP MASS OF IMP ACTOR 

C0NM2. 200. 999. 0.9. 587-5. 0.0. 0.0. 00. . +CON2-2 
4C0N2— 2.3.745-6. .3.745-6, . .3.745-6 
$ MATERIAL PROPERTIES 

PROD .1,11.0.196.6.1 4-3 .0.25 
MAT1 .11.10.0+6, .0.33.2.5-4, . . .+MAT1-1 
+MAT1-1 , 35 . 0E6 .36 . 0E6 , 27 .0E6 
$ BOUNDRY CONDITIONS 

SPC1 . 4 . 23456 . 1 . THRU , 1 45 

$ REMOVE DEGREES OF FREEDOM FROM IMPACTOR 

SPC1 , 4. 23456. 999 

$ TRANSFER FUNCTION TO DEFINE INDENTATION 

TF.111 .1001 .0.+1 .0,0. 0,0.0. . ,+TF-l 

+TF-1 .999.1.-1 .0.0. 0.0.0. . . ,+TF— 2 

+TF— 2 .1,1, 1.0, 0.0, 0.0 

$ TIMING 

TSTEP . 7 . 2500 . 2 . 0-7 . 25 

s LOAD DEPENDENT ON DISPLACEMENT OF IMPACTOR 

N0LIN3. 5. 1. 1. 6.24+6. 1001. 1. 1.5 
$ SLOW DOWN IMPACTOR 

N0L1N3. 5. 999. 1. -6.24+6. 1001. 1. 1.5 

$ INITIAL CONDITIONS: IMPACTOR VELOCITY - 62.1 IN/SEC 

TIC. 1 .999.1 .0.0.62.1 

ENDDATA 
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ANA LEX - STRUCTRAL MECHANICS 8RANCH 


ID IM»ACT. PLATE 
APP OISP 
TIME 120 
SOL 27 
CEND 

TITLE - IMPACT OF PLATE 49X49 : CENTERED ELEMENT 
SUBTITLE - TRANSIENT ANALYSIS: FIXED-FIXED: NO SYMMETRY 
SPC ■ 1 
IC - 3 

NONLINEAR - 5 
TSTEP - 1 
TFL - 111 

SET 15 - 999 . 2525,2526,2625,2626 
NLLOAO - 15 

SET 20 - 2525. 3325. -41 25 
STRESS - 20 
BEGIN BULK 

$ .... EXTRA POINT TO HOLD INDENTATION ................. 

EPOINT. 10001 

S «... IMP ACTOR 3/8 IN DIAMETER 

GRID. 999, .0.0.0. 0,-0. 1875 

CONM2. 200. 999. 0.8. 096-5. 0.0. 0.0, 0.0, , +CON2-2 
4CON2-2 . 7 . 459-6 . . 7 . 459-6 . . . 1 . 423-6 


$ ...... GRIDS AND CQUAD4 ELEMENTS DEFINING THE PLATE GO HERE 

$ MATERIAL PROPERTIES... MAT2 CARDS GENERATED BY C06STRAN 


PSHELL. 1.101. 0.15. 201, 1.0 

MAT2. 101 . 4. 3E+06 .2 .9E+05 1 . 7E-03 . 2 .8E+06 ,-3 . 4E-02 , 5 . 7E+05. 1 8E-04 +A101 
+A101 .5.8E— 06.8. 9E— 06,5. 0E— 13 

MAT2.201 .5. 7E+06 . 2 . 9E+05,— 1 . 9E— 04 . 1 . 4E+06 .-3 .8E-03 . 5.7E+05 
$ BOUNDRY CONDITIONS 

SPC1. 1. 123456. 101, THRU. 150 

SPC1 . 1. 123456. 5001. THRU. 505e 

SPC1 , 1. 123456. 101 

-.-.-.*100 

-48 

SPC1 , 1. 123456, 150 

-.-.-.•100 

-48 

SPC1 , 1. 12456. 999 


GRDSET 6 

$ TIME STEP INFO 

TSTEP. 1.2000. 1.0-7. 10 

S LOAD DEPENDENT ON RELATIVE DISPLACEMENT OF IMPACTOR 


NOLIN3. 5.2525, 3. +1.945+5, 10001, 0. 1.5 

NOLIN3. 5.2526. 3. +1.945+5, 10001, 0. 1.5 

NOLIN3. 5.2625. 3. +1.945+5, 10001. 0. 1.5 

NOLIN3. 5,2626, 3. +1.945+5. 10001. 0. 1.5 

J SLOW DOWN IMPACTOR 

NOLIN3. 5. 999. 3,-7.779+5, 10001, 0, 1.5 

S TRANSFER FUNCTION TO CALCULATE INDENTATION 

TF.111 .10001 .0.+1 .0,00.0.00.0, . ,+TF— 1 

+TF-1 .999.3.-1 .0.00.0.00.0, , . ,+TF— 2 

+TF-2. 2525. 3. +0.25, 00. 0.00.0, . . .+TF-3 

+TF-3. 2526. 3. +0.25, 00. 0,00.0. , , .+TF-4 

+TF-4. 2625. 3. +0.25, 00. 0,00.0, . . .+TF-5 

+TF-5. 2626, 3. +0.25. 00. 0,00.0 

S INITIAL CONDITIONS: IMPACTOR VELOCITY - 1470 IN/SEC (122.5 FT/SEC) 

TIC. 3. 999. 3. 0.0. 1470.0 J 

ENOOATA 


134 



ABSTRACT 

The propagation of power flow through a dynamically 
loaded beam model with 90 degree bends is investigated using 
NASTRAN and McROW. The transitioning of power flow 
types (axial, torsional, and flexural) is observed throughout 
the structure. To get accurate calculations of the torsional 
response of beams using NASTRAN, torsional inertia effects 
had to be added to the mass matrix calculation section of the 
program. Also, mass effects were included in the calculation 
of BAR forces to improve the continuity of power flow 
between elements. The importance of including all types of 
power flow in an analysis, rather than only flexural power, is 
indicated by the example. Trying to interpret power flow 
results that only consider flexural components in even a 
moderately complex problem will result in incorrect 
conclusions concerning the total power flow field. 

INTRODUCTION 

Methods for calculating power Hows in dynamically loaded finite 
element models using NAS IRAN (Rigid Format 8 - Direct Frequency 
Response) and McPOW (Mechanical POWer) were developed previously. 1 
The power flow equations for beam elements derived in that paper included all 
forms of dynamic energy propagation: flexural, longitudinal (or axial), and 
torsional. The flexural waves were split into shear and moment components. 

The majority of procedures employed in other studies (see the list ol 
references in Hambric 1 ) only consider flexural vibration in their calculations ol 
power flow. This can be dangerous if an analyst is investigating the energy 
propagation characteristics of a complex structure. Though flexural vibration 
is in most cases the dominant response in a dynamically excited beam, 
different kinds of propagation will occur in structures with even a small degree 
of complexity, such as a simple beam model with 90-degree bends. 
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Such a model is tested here using a frequency range spanning several 
resonances and types of motion. Plots showing the contributions of the 
different forms of power How to the total power travelling through the system 
are shown, and illustrate the importance of all types of energy propagation to 
the power flow method. 

To improve the accuracy of both the finite element solution and the 
power flow solution ol the problem, a few modifications were made to 
NASTRAN and MePOW. First, to show the importance of torsional power 
flow, a capability to calculate dynamic torsional forces and corresponding 
angular velocities is required. Therefore, torsional inertias were added to the 
coupled mass matrix formulation of the BAR element. Also, since the beam 
element force calculation algorithm in NASTRAN considers only stiffness 
effects, mass and damping effects had to be added to MePOW to modify the 
element forces. 

MKTHODOLOGY 

The procedure for solving for the power flow field in a finite element 
model using NASTRAN and MePOW is: 

1. Run Rigid Formal 8 (Direct Frequency Response) on a NASTRAN data 
deck (using the ALT HR statements shown in Ref. 1 to output force and 
velocity data blocks to the OUTPUT2 file). Coupled mass formulations 
should always be used. 

2. Run MePOW using the binary data in the OUTPUT2 file as input. 

General Methods 

A typical power flow cycle is shown in Fig. 1. The figure shows an 
arbitrary structure mounted to a connecting structure by a spring and clamper 
coupling. A dynamic load is applied, and energy flows into the structure at the 
load point. 'Pile input power then flows through the structure along multiple 
flow paths denoted by arrows whose lengths represent power flow magnitudes. 
As the energy flows toward the mounting, it is dissipated by material damping 
and sound radiation into a surrounding medium, and the flow arrows shorten. 
The flow and dissipation processes continue until the remaining energy exits 
the structure through the mounting and flows into the connecting structure. 
Though only one power entry point and one exit point are shown in this 
drawing, multiple loads and mountings may exist. A classic text which 
describes the flow of structure-borne sound is the book by Cremer, Ileckl, and 
IJngar. 2 

The structural dynamics problem may be solved using NASTRAN. The 
structure may be modeled with various element types; mountings are modeled 
with scalar spring, damping, and mass elements. Constraints and loads are 
directly applied. The steady-stale response for the model is solved for a given 
excitation frequency, and the power flow variables are calculated. 
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Fig. 1. Sample Power Flow Diagram. 


Power is defined as the time-averaged product of a lorce with the in- 
phase component of velocity in the direction of the force. I ; or time-harmonic 
analysis, where complex numbers are used, this calculation may be visualized 
as taking the dot product of the force and velocity phasors. (There is no factor 
1/2 in the following power equations il the assumption that lorces and 
velocities are "effective , ‘ values rather than amplitudes is made. With this 
assumption, consistency is maintained, and there is no mixing ol effective and 
peak quantities in this formulation.) 

Multiplying one complex number by the in-phase part ol another 
complex number is the same operation as multiplying the first number by the 
complex conjugate of the other number and taking the real part ol the result. 
Therefore a general formula for power How in a structure is 

Power = Real [Fv |, (0 

where 

F = force, and 

v* — complex conjugate of velocity. 


137 



Power Flow Equations 


The equations for power flows in BAR elements are repeated here. A 
diagram of the BAR element and its NASTRAN force output conventions is 
shown in Fig. 2, where Plane 1 is vertical and Plane 2 is horizontal. 



Fig. 2. The BAR Element 


Since a beam is a one-dimensional element, energy Hows in only one direction: 
in the local x direction, or along the length of the beam. The total power How 
for a beam element is 

P x = Real [ - ( I v x v x T V i v y + V 2 v z + 1 — M 2 0 * -f M j 0* ) ] , (2) 

where 

F x = axial force, 

V ! = shear force in y direction, 

V 2 = shear force in z direction, 

T = torsion about x, 

= bending moment about y, 

M] = bending moment about z, 

V; = translational velocities in direction i, and 
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S\ = rotational velocities about axis i. 

The negative sign on the result comes from force and displacement direction 
conventions for the element. The negative sign on the M 2 term reflects the 
NASTRAN force output convention. In Fig. 2, M 2 is shown as positive in the 
opposite sense to 6 y . Therefore, M 2 0* is opposite in sign to the other power 
llow components. 

NASTRAN Modifications 

Torsional Inertias 

NASTRAN currently does not consider torsional inertias in its beam 
element formulation. Therefore, all torsional results (angular displacements 
and torques) arc based on stiffness only, and are essentially those of a static 
problem solution. To remedy this, torsional inertias were added to the 
coupled mass formulation. At the point in NASTRAN where the basic 
element mass matrix is formed, no consideration is given to beam offsets or 
beam orientation; all mass coefficients (as well as stiffness) are calculated in 
the local beam coordinate system. 

The torsional mass moment of inertia of a beam is p L J x /2, where p is 
the mass density, L is the beam length, and J x is the polar area moment of 
inertia. In the standard consistent mass matrix for a beam/ this value is 
broken up into 2/3 and 1/3 components; 2/3 of the value is placed at the 
diagonal, and 1/3 is placed at the coupled degree of freedom (the node on the 
other end of the beam). The same fractions are used for the translational, or 
axial masses. In NASTRAN, however, the coupled mass formulation uses an 
average of lumped and consistent formulations to reduce error. This average 
changes the components to 5/6 and 1/6 of the total value. Since these values 
are currently used for the axial masses in NASTRAN, they were also used for 
the torsional inertias. 

Element Force Calculations 

NASTRAN element forces are currently calculated by multiplying 
element stiffness matrices by element displacement vectors. Both damping 
and mass effects are ignored. The damping in a stiffness element is actually in 
the form of a loss factor, which generates a complex stiffness matrix. All 
stiffness terms are multiplied by 1.0 + if/. For most dynamic analyses, 
neglecting the b] term is acceptable since it is generally small. For a power 
flow analysis of a highly reverberant structure, however, ignoring the loss 
factor is disastrous. In a highly reverberant structure, the force and velocity at 
a given point are close to 90 degrees out of phase. Since power llow is defined 
as the dot product of these two components, a small change in the phase of 
the force has large effects on the calculated element powers. 

Neglecting the element mass matrices, whose components are several 
orders of magnitude less than those of the stiffness matrices, has less drastic 
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effects on the power How solution, since at low frequencies the masses will 
have little effect on the force calculations (the element mass matrix is 
multiplied by — c o 2 to take the second time derivative of the corresponding 
displacements). However, when high frequency analyses are performed on a 
model, the — e o 2 multiplying factor becomes more significant, and neglecting 
the mass contributions will cause some error in the force calculations. Errors 
in element forces cause errors in element power Hows. 

Including these missing effects in NASI RAN is complicated by the fact 
that the element force calculation algorithm splits the problem into real and 
imaginary parts. The element stiffness matrices are multiplied by the real 
parts of the displacement vectors to calculate real force components, and the 
process is repeated for the imaginary components. Adding an imaginary term 
to the stiffness matrices causes new terms to be generated in the multiplication 
(imaginary stiffness x imaginary displacement and imaginary stiffness x real 
displacement). There is also no frequency dependence in the current 
algorithm, since stiffness are frequency independent. Mass matrices, however, 
must be multiplied by the — ca 2 term mentioned above, so they must be 
recalculated for every frequency. 

To avoid these complications, the element force calculations were 
temporarily moved to McPOW. The element mass and complex stiffness 
matrices are recalculated on a local element level, and combined with local 
element displacements to solve for element forces. A force vector with 12 
entries is the result; shears in the local y and z directions, moments about the 
local y and z directions, axial forces, and torques are solved for at each grid 
point. In NASTRAN only eight forces are calculated, because only moments 
are calculated at both ends of a beam element. Beam power Hows are 
therefore calculated at each end of the element using only the forces at that 
end and the corresponding grid velocities. The average of the powers at the 
ends is taken to find an element power tlow. 

TKS1 PROBLEM 


Problem Statement 

The beam model that was analyzed is shown in Fig. 3. All three 
sections have the same cross section and material properties. Dashpots 
(DAMP2 elements) of value \0° were applied at the model’s end in all six 
degrees of freedom. A unit load was applied at the top end of the model in 
the longitudinal direction (along the - z axis) over a frequency range of 1 to 250 
Hz swept in 1 Hz increments. The finite element model consists of 152 
elements and 153 grid points. Grid and clement numbering starts at the left 
end of Link 1 and proceeds up to the end of Link 3. 

Using the local beam element coordinate systems shown in lug. 3, the 
following table of force balances at the corners (Lank 3 to Link 2, Link 2 to 
Link 1) was generated. 
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Fig. 3. Test Problem Geometry 


Link 3 

Link 2 
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The subscripts on the shears and moments refer to the plane in which the 
forces occur (sec Fig. 2). This table can be used to track the propagation ot 
power How through the structure. For example, the longitudinal power input 
to Link 3 will travel down the beam in axial waves to the first bend and 
become shear power How in the z direction in Link 2. This shear power will 
interchange with moment power along the beam (the sum of the shear and 
moment components is the total flexural power flow in the beam). Any shear 
power that exists at the lower end of Link 2 will transition to more shear 
power in Link 1 . 


141 







75.0 200.0 225.0 





Results and Discussion 

The computed power input curve over the excitation frequency range is 
shown in Fig. 4. The power input peaks correspond to various resonances in 
the structure. Most are flexural, but some axial and torsional inodes influence 
the power input curve. The longitudinal modes of Link 3 cause power input 
peaks (at 190 Hz and above), as well as the torsional modes of lank 1 (at 151 
Hz and above). 

In this model, the power flow path is independent of frequency. The 
total power must always flow from the input point at the end oi Link 3 to the 
dampers at the beginning of Link 1. This simplifies the interpretation of the 
results, since the directions of total power flow are established. 

The types of power flow in a given link are not so well-defined. 
Whether the dominant path in a link is flexural, axial, or torsional, depends on 
the motion of the structure. Fig. 5 shows the two most common types of 
motion paths for this problem. The displacement field of Diagram 1 occurs 
most often. The axial load applied to Link 3 drives the entire structure 
forward and backward over a frequency cycle. The dominant power flow in 
Link 3 is axial; the dominant power flow in Link 2 is flexural; and torsional 
and flexural power flows are dominant in Link 1, since the input load applies 
both a torque and bending moment to the link. 

In Diagram 2 of big. 5 a different type of motion is shown. The axial 
load still drives the upper half of the structure in the same direction, but the 
lower half moves in the opposite direction. This type of motion is not what 
one would expect in a static problem, blit the dynamic characteristics ol the 
structure produce this type of motion in various frequency ranges. 

Due to this motion path, the axial power flow travelling down Link 3 
becomes flexural, torsional, and axial in Link 2. The torsional and axial 
components appear because the link is twisted and stretched by the opposite 
directions of motion of the two ends. The torsional power in Link 2 becomes 
flexural power in Plane 1 in Link 1, and the axial power in Link 2 turns into 
flexural power in Plane 2 in Link 1. The flexural power in Link 2 becomes 
torsional and flexural power in Link 1 as before (Diagram 1). 

Considering these modes of power transitioning, the power flow plots in 
Figs. 6-8 may be interpreted. Each plot shows the contributions of flexural, 
torsional, and axial power flow as a percentage of the total power flow in the 
center of each link. 

Fig. 8 shows the power components in Link 3. Since the input power is 
in the longitudinal direction, the majority of the power in this link is axial. At 
certain frequencies, the percentage of axial power is greater than 100 percent. 
The large axial percentage arises because at certain frequencies, reflected 
waves carry power in the opposite direction (toward the load). Three flexural 
resonances in the structure cause reflected power just before 50 Hz, along with 
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Fig. 5. Dominant motion paths for test problem 


five others right after 100 Hz. Between 200 and 250 Hz, some flexural and 
torsional resonances cause more reflected powers. 

Fig. 7 shows the power components in Link 2. The dominant type of 
power is the flexural component in Plane 1, and is denoted by the solid curve. 
This type of power field corresponds to the motion type shown in Diagram 1 in 
Fig. 5. However at certain frequencies, the power flow pattern of Diagram 2 
becomes dominant, and axial and torsional power become important. In most 
cases, the axial power flows forward (away from the load point), and the 
torsional power is backward (reflected toward the load point). These 
tendencies occur at the same frequencies as the reflected power waves do in 
Link 3 (shown in Fig. 8). This behavior indicates that the flexural power in 
Plane 1 and the torsional power cause reflected flexural powers in Planes 2 and 
1 respectively in Link 3. 

Fig. 6 shows the power distribution in Link 1. In this case, all power 
components are positive, implying that the reflected power waves in Links 2 
and 3 originate from the joint connecting Links 1 and 2. In Link 1, flexure in 
Plane 1 and torsion are the dominant components of power flow. Flexural 
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Power How types in Link 1. 
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Power flow types in Link 3. 



motion in Plane 2 and axial motion cause power peaks at the same frequencies 
observed in I igs. 7 and 8, indicating the type of motion shown in Diagram 2 
in 1 i.,. 5. A toisional mode in Dink ] accounts for the peak in the torsion 
curve at 150 Hz, along with an input power peak at the same frequency (see 
Pig. 4). 

In spite ot the large variation in percentages of power types in the plots, 
all the power curves add up to 100 percent, as expected. In addition, the total 
power llow in the structure at all frequencies is at a maximum at the load 
point, and smoothly decreases to a minimum at the connection point to the 
dampers. I he steady decrease in power is due to structural damping. The 
remaining power is all dissipated by the connected dampers. 

I his example illustrates the importance of all types of power 
components in a power llow analysis. Imagine trying to discern a meaningful 
power llow field from only flexural powers in this example. The detected 
powers in l ink 3, which is adjacent to (he input load, are all in the oposite 
direction, or toward the load. In Link 2, the analyst would see a sudden jump 
m power to values that arc higher than that of the input power. Finally, in 
Link 1, sporadic power curves with values near the input power at frequencies 
below 100 IIz and values near zero alter 100 Hz would be found. Confusion 
would suiely result, with erroneous conclusions soon following. Difficulties 
like these would be compounded in a real application with some degree of 
complexity. 


CONCLUSIONS 

The modifications made to NASTRAN and MePOW are critical to the 
power llow method. Without torsional inertias applied to the beam element 
mass matrices, any torsional effects in a dynamic problem are static. None of 
the torsional power flows present in the example problem would exist, causing 
incorrect total power flow fields. Adding mass and damping effects to the 
element loree calculation algorithm is also important. In a reverberant 
structure where forces and velocities are nearly 90 degrees out of phase with 
each other, accurate calculations are necessary to guarantee good power flow 
results. A small change in the phase of an element force, caused by neglecting 
the material loss factor, could cause large errors in element power" flows & 
Also, at higher frequencies, element mass terms can become significant and 
nllect the element force magnitudes, and hence the element power magnitudes. 

The addition of torsional inertia to the beam element mass matrix 
lormulation was straight-forward. The addition of damping and mass effects to 
the element force calculation routines, however, was almost impossible. In 
fact, the changes had to be made to MePOW instead of NASTRAN. The 
implementation difficulties were due to the way NASTRAN handles complex 
analysis: the solutions are broken into real and imaginary parts. When the 
pioginm was in its formative stages, UNI VAC' computers were supported. 
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The UNIVAC, unfortunately, had no way of handling double precision 
complex arithmetic. Therefore, no complex numbers or FORTRAN complex 
functions are used in the element force calculation sections of the program. 
With this approach, a simple complex calculation like 
[—or [M] e + (1 +i//)[K] e j {d} e must be split up into four calculations. Also, 
since the calculation is frequency-dependent, the NASTRAN element force 
subroutines are not currently able to handle it. Since the UNIVAC has all but 
disappeared from the COSMIC NASTRAN computing arena and most 
modern computers support double precision complex arithmetic, perhaps the 
way NASTRAN handles complex problems should be modified. 

The importance of including longitudinal and torsional components with 
flexural ones in a power flow analysis was shown in the example problem. 
Measuring flexural power alone will not give an accurate indication of the total 
power flow field in even a marginally complex problem. In the case of the 
example problem, reflected flexural waves actually indicated a reversal of 
power flow in the model, where the direction of flexural power was in the 
opposite direction of the input power. Trying to interpret power flow results 
that only consider flexural components will result in incorrect conclusions 
concerning the total power llow field. 
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ABSTRACT 

NASHUA is a coupled finite element/boundary element capability built around 
NASTRAN for calculating the low frequency far-field acoustic pressure field radiated or 
scattered by an arbitrary, submerged, three-dimensional, elastic structure subjected to either 
internal time-harmonic mechanical loads or external time-harmonic incident loadings. This 
paper describes the formulation and use of NASHUA for solving such structural acoustics 
problems when the structure is fluid-filled. NASTRAN is used to generate the structural 
finite element model and to perform most of the required matrix operations. Both fluid 
domains are modeled using the boundary element capability in NASHUA, whose matrix 
formulation (and the associated NASTRAN DMAP) for evacuated structures can be used 
with suitable interpretation of the matrix definitions. After computing surface pressures and 
normal velocities, far-field pressures are evaluated using an asymptotic form of the 
Helmholtz exterior integral equation. The proposed numerical approach is validated by 
comparing the acoustic field scattered from a submerged fluid-filled spherical thin shell to 
that obtained with a series solution, which is also derived in this paper. 


INTRODUCTION 

Two basic problems in computational structural acoustics are (1) the calculation of the 
acoustic pressure field radiated by a general submerged three-dimensional elastic structure 
subjected to internal time-harmonic loads, and (2) the calculation of the acoustic pressure 
scattered by such a structure subjected to an incident time-harmonic wavetrain. The most 
common, as well as the most accurate, approach for solving these problems at low 
frequencies is to couple a finite element model of the structure with a boundary element 
model of the surrounding fluid. This is the approach taken by NASHUA, which is a 
boundary clement program built around NASTRAN, a widely-used finite element computer 
program for structural dynamics. 

Several previous papers (Ref. 1-4) described the basic formulation and development for 
acoustic radiation and scattering from evacuated structures. Here we describe the 
lormulation and use of NASHUA for modeling submerged structures which are fluid-filled. 
Internal fluid can occur because the structure is free-flooded or contains fluid-filled tanks. It 
is possible to use existing NASTRAN capability to model the interior fluid with finite 
elements (Ref. 5-7), but three-dimensional models with large numbers of fluid degrees of 
freedom might result. An attractive alternative to the fluid finite element model is to 
represent the contained fluid using a boundary element approach. In principle, any computer 
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program capable of generating the appropriate boundary element matrices for an exterior 
fluid is also capable of generating such matrices for the complementary region (the interior 
region). NASTRAN’s versatility in user-controlled matrix operations (DMAP) makes the 
implementation of such an approach straightforward. 


THEORETICAL APPROACH 

The basic theoretical development for NASHUA’S radiation and scattering approach 
for evacuated structures has been presented in detail previously (Ref. 1-4). For 
completeness, this paper summarizes that approach and describes the changes necessary to 
model the interior fluid with boundary elements in the same procedure. There is no 
requirement that the interior and exterior fluids be the same. 


The Surface Solution for Evacuated Structures 

Consider any submerged three-dimensional, evacuated elastic structure subjected to 
either internal time-harmonic loads or an external time-harmonic incident wavetrain. If the 
structure is modeled with finite elements using NASI RAN, the resulting matrix equation of 
motion can be written as 

Zv = F -GAp, W 

where matrix Z (of dimension s x s) is the structural impedance, vector v (s x r) is the 
complex velocity amplitude for all structural DOF (wet and dry) using the coordinate systems 
selected by the user, vector F (s x r) is the complex amplitude ol the mechanical forces 
applied to the structure, matrix G (s x f) is the rectangular transformation of direction 
cosines to transform a vector of outward normal forces at the wet points to a vector of 
forces at all points in the coordinate systems selected by the user, matrix A (f x f) is the 
diagonal area matrix for the wet surface, and vector p (f x r) is the complex amplitude of 
total pressures (incident + scattered) applied at the wet grid points. In this equation, the 
time dependence exp(iwt) has been suppressed. In the above dimensions, s denotes the total 
number of independent structural DOF (wet and dry), f denotes the number of fluid DOb 
(wet points), and r denotes the number of load cases. In general, the surface areas, the 
normals, and the transformation matrix G are obtained in NASHUA from the NASTRAN 
calculation of the load vector resulting from an outwardly directly static unit pressure load on 
the structure’s wet surface. 

In Eq. 1, the structural impedance matrix Z, which converts velocity to force, is given 

by a 

Z = (-arM + iwB + K)/(ice), ( 2 ) 

where M, B, and K are the structural mass, viscous damping, and stiltncss matrices, 
respectively, and u > is the circular frequency of excitation. For structures with a nonzero loss 
factor K is complex. In addition, K can include the differential stiffness effects of 
hydrostatic pressure, if any (Ref. 3). A standard NASTRAN finite element model of the 
structure supplies the matrices K, M, and B. 

For the exterior fluid domain, the total fluid pressure p satisfies the Helmholtz, 
differential equation 
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v 2 p + k 2 p = 0, (3) 

where k = cu/c is the acoustic wave number, and c is the fluid sound speed. Equivalently, for 
the exterior fluid, p is the solution of the Helmholtz integral equations (Ref. 8) 


4 - 4 c i( x ) D ( r )Js 


p(x')/2 - pi, 
• P(x') - pi, 
-Pi> 


x' on S, 
x' in E, 
x' in I, 


( 4 ) 


where S, E, and I denote the surface, exterior, and interior domains, respectively, p T is the 
incident free-held pressure (if any), r is the distance from x to x' (Fig. 1), D is the free-space 
Green’s function 


D(r) 



q 


dp_ 

dn 


-iw/>v n , 
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Fig. 1. Notation for Helmholtz Integral Equation 


p is the fluid mass density, and v n is the outward normal velocity on S. As shown in Fig. 1, x 
in Eq. 4 is the position vector for a typical point Pj on the surface S, x' is the position vector 
for the point Pj on the surface or in the exterior field, the vector r = x' - x, and n is the unit 
outward normal at Pj. We denote tire lengths of the vectors x, x', and r by x, x’, and r, 
respectively. The normal derivative of the Green’s function D is (Ref. 1) 


dD(r) 

<9n 


(7) 


where 8 is the angle between the normal n and the vector r, as shown in Fig. 1. 

All three integral equations in Eq. 4 are needed for exterior fluids. The surface 
equation provides the fluid impedance at the fluid-structure interface. Since the surface 
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equation exhibits non-uniqueness at certain discrete characteristic frequencies (Ref. 9), the 
interior equation is used to provide additional constraint equations which ensure the required 
uniqueness. The exterior equation is used to compute the exterior pressure field once the 
surface solution (which includes the fluid pressure and its gradient) is known. 


The substitution of Eqs. 6 and 7 into the surface equation (4) yields 

- f p(x) — A (ik + — ) cos p dS = i up f v„(x) ^ dS + Pi, 

2 J s 1 v ' 47iT r 4;Tr 


x' on S. (8) 


This integral equation relates the total pressure p and normal velocity v n on S. If the 
integrals in Eq . 8 are discretized for numerical computation (Ref. 1), we obtain the matrix 
equation (for the exterior fluid) 

Ep = Cv„ + pi , ( 9 ) 

where vector p (of dimension 1 x r) is the vector of complex amplitudes ol the total pressure 
on the structure’s wet surface, matrices E and C (both 1 x f) arc lully-populated, complex, 
nonsymmetrie, and Irequency-dcpendent, and vector pj (f x r) is the complex amplitude ol 
the incident pressure vector. The number ol unknowns in this system is 1, the number ol 
wet points on the fluid-structure interface. 

The normal velocities v n in Eq. 9 are related to the total velocities v by the same 
rectangular transformation matrix G: 

v n = G t v, ( 10 > 

where T denotes the matrix transpose. If velocities v and v„ are eliminated from Eqs. 1, 9, 
and 10, the resulting equation for the coupled fluid-structure system is 


(E + CG r Z -1 GA) p = cg t z- 1 f + p,. (ID 

This equation is solved for the total surface pressures p, since the rest of the equation 
depends only on the geometry, the material properties, and the frequency. Since the two 
right-hand side terms in Eq. 11 correspond to mechanical and incident loadings, only one ol 
the two terms would ordinarily be present for a given case. The details ol the incident 
pressure vector pj for scattering problems were presented previously (Ret. 2). and will not be 
repeated here. 

The velocity vector v for all structural DOF is recovered by solving Eq. 1 for v: 

v = Z _1 E -Z _1 GAp. (12) 

The surface normal velocity vector v n is recovered by substituting this solution lor v into Eq. 

10 . 


Modeling Interior Fluid 

The theoretical development presented in the preceding section can be modified slightly 
to account also for an interior fluid. The wave equation, Eq. 3, applies also to interior fluids. 
Although all three integral equations in Eq. 4 are generally needed for exterior fluids, only 
the surface equation is needed to represent the surface impedance ol interior fluids. Eq. 4a 
also applies to interior fluids if the incident pressure p] is set to zero, and the normal vector 
n is negated. That is, the surface integral equation applies to both exterior and interior fluids 
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so long as the unit normal is always directed from the structure into the fluid. One other 
consideration, perhaps unique to NASHUA, is that wet surface curvatures (which are 
needed in the calculation of the "self” terms in matrix E) are negative at some interior points 
(Ref. 1). 

A matrix equation similar to Eq. (9) is therefore obtained for the interior fluid except 
that the incident pressure p x is zero. The fluid matrices E and C are different for exterior 
and interior domains (even if the separating surface S has infinitesimal thickness) because 
the normals are of opposite sign. 


Two-Fluid Formulation 


Denote the exterior fluid as Fluid 1 and the interior fluid as Fluid 2, and use the 
subscripts 1 and 2 to refer to these two domains. Also define new pressure and normal 
velocity unknowns p and v n which include the solutions for both fluid domains: 


P = 




( 13 ) 


Since there is no direct fluid coupling between the interior and exterior fluids, and the 
incident pressure vanishes in the interior domain, Eqs. 1, 9, 10, and 11 apply also to the 
two-fluid problem if the new definitions in Eq. 13 are used, and the matrices A, G, E, C, 
and pi are re-defined as 


A = 



G 


G, G 2 , 




(14) 


The principle benefit of formulating the two-fluid problem in this way is that the required 
modifications to extend the procedure to three or more independent fluid domains is now 
clear. 


The Far-Field Calculation 

With the solution for the total pressures and velocities on the surface, the exterior 
Helmholtz integral equation, Eq. 4b, can be integrated to obtain the radiated (or scattered) 
pressure at any desired location x' in the exterior field. We first substitute Eqs. 5-7 into Eq. 
4b to obtain 

1 e -ikr 

p( x ') = I f iw /> v n( x ) + ( ik + — )p( x ) cos p\ ~~ A dS > x' ill E, (15) 

s r 4 TIT 

In applications, however, the field pressures generally of interest are in the far-field, so we 
use instead the asymptotic form of Eq. 15 (Ref. 1): 

p(x') = — j — — I [/>cv n ( x ) + p(x) cos /?]e’ kx cos “ dS, x' in E, x' » d, (16) 

47TX S 

where d is a characteristic dimension of the structure, and a is the angle between the vectors 
x and x' (Fig. 1). For far-field points, cos ft is computed using the asymptotic approximation 
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For both Eqs. 15 and 16, numerical quadrature is used. 


OVERVIEW OF SOLUTION PROCEDURE 

The NASHUA solution procedure uses NASTRAN to generate the matrices K, M, B, 
and F and to generate sufficient geometry information so that the matrices E, C, G, A, and 
Pi can be computed by a separate program called SURF. Then, NASTRAN DMA? is used 
to form the matrices appearing in Eq. 11, which is solved for the total pressures p (in both 
fluid domains) using the block solver OCSOLV (Ref. 10). Next, NASTRAN DMAP is used 
to recover the surface normal velocities v n and the vector v of velocities at all structural 
DOF (NASTRAN’s "g-set"). This step completes the surface solution. Then, with the total 
pressures and velocities on the (exterior) surface, the asymptotic (far-ficld) form of the 
Flelmholtz exterior integral equation is integrated in program FAROUT to compute the far- 
field radiated pressures. Various tables and graphical displays are generated. 

The overall setup of the solution procedure is organized into four steps. In Step 1, a 
separate NASTRAN structural model is prepared and run for each unique set of symmetry 
constraints and each fluid region. Since, for general three-dimensional analysis, up to three 
planes of reflective symmetry are allowed, there would be one, two, four, or eight such runs 
for each fluid region. Since the purpose of this step is to generate a file containing geometry 
information and a checkpoint file for subsequent use in the other steps, the only difference 
between the two runs associated with a given symmetry case is the specification of the 
outwardly directed unit pressure load which defines the wet surface for a given fluid region. 

For each symmetry case and drive frequency, several programs are run sequentially to 
form Step 2. For each fluid region, the SURF program reads the geometry file generated by 
NASTRAN in Step 1 and, using the Helmholtz surface and interior integral equations, 
generates the fluid matrices E}, E 2 , C[, and C 2 , the area matrices A\ and A 2 , the structure- 
fluid transformation matrices Gi and G 2 , the incident pressure vector pn, and a geometry 
file to be used later by the far-field integration program FAROUT in Step 3. In addition, a 
partitioning vector is generated to facilitate the merging and partitioning of the various 
matrices associated with the two fluid domains. 

The two SURF jobs in Step 2 are followed by a NASTRAN job which takes the 
structural matrices K, M, B, and F from Step 1 and the matrices generated by the SURF 
jobs and forms the matrices in Eq. 11, where the definitions in Eq. 14 apply. Eq. 11 is then 
solved for the total surface pressure vector p by program OCSOLV, which is a general out- 
of-core block solver designed specifically for large, full, complex, nonsymmetric systems of 
linear, algebraic equations. NASTRAN is then re-entered in Step 2 with p so that the 
velocities v and v n can be recovered using DMAP operations. The surface pressures, 
normal velocities, and full g-set displacements are then reformatted, sorted, and merged into 
a single file (for each symmetry case) using program MERGE. Recall that there are one, 
two, four, or eight possible symmetry cases. 

Steps 1 and 2 arc repeated for each symmetry case. After all symmetry cases have 
been completed and merged, program FAROU T (Step 3) combines the symmetry cases and 
integrates over the surface. The far-field pressure solution is obtained by integrating the 



surface pressures and velocities using the asymptotic (far-field) form of the exterior 
Helmholtz integral equation, Eq. 16. Output from FAROUT consists of both tables and files 
suitable for various types of plotting. 

The remaining steps in the NASHUA procedure are for graphical display. Deformed 
structural plots of the frequency response are obtained by restarting NASTRAN (Step 4) 
with the checkpoint file from Step 1 and a results file from FAROUT. In addition, animated 
plots can be generated on the Evans & Sutherland PS-330 graphics terminal using the 
CANDI program written for the DEC/VAX computer by R.R. Lipman of DTRC (Ref. 11). 
X-Y plots of various quantities (both surface and far-field) versus frequency may be obtained 
using IPLOT or other interactive plotting programs (Ref. 12). Polar plots of the far-field 
sound pressure levels in each of the three principal coordinate planes can also be generated 
using the interactive graphics program FAFPLOT (Ref. 1). 


DMAP ALTERS 

Several DMAP alters are used in the overall NASHUA procedure to implement the 
preeedure described in preceding section. For Step 1, the following alter is used: 


ALTER 

ALTER 

ALTER 

GP3 

ALTER 

SSG1 

SSG2 

OUTPUT2 

OUTPUT2 

OUTPUT2 

PARAMR 

COND 

PARAMR 

DIAGONAL 

ADD 

RBMG2 

SSG3 


SDR1 

TA1 

DSMG1 

EQUIV 

COND 

MCE2 

LABEL 

EQUIV 

COND 


1 $ NASHUA STEP 1, COSMIC 1988 RF8 (REVISED 12/7/89) 

2,2 $ DELETE PRECHK 
21,21 S REPLACE GP3 

GEOM3,EQEXIN,GEOM2/SLT,GITT/S,N,NOGRAV/NEVER=l $ 
117,117 S REPLACE FRRD 

SLT.BGPD r,CSTM,SIL,EST,MPT,GPTT,EDT,MGG,CASECC,DIT,/ 
PG , , ,,/LUSET/N SKIP $ PG 

USET,GM,YS,KFS,GO,DM,PG/QR,PO,PS,PL S PL 
AXIC,BGPDT,EQEXIN,USET,PG $ 

PL.CSTM ,ECT, , $ 

,,„//- 9 $ 

//*EQ*//C,Y,IISP=0./0.////NOHSP S 
LBL4D,NOIISP S SKIP DIFF. STIFF. IF NO HYDRO. P 
//*COMPLEX//C,Y,HSP=0./0./HSPC $ HSP+I*0 
KAA/KDIAG/*SQUARE*/1.0 $ 
KAA,KDIAG/KAAD/(1.0,0.0)/(1.E-6,0.) S 
KAAD/LLL S FACTOR KAA 

LLL,KAAD,PL,LOO,KOO,PO/ULV,UOOV,RULV,RUOV/OMIT/ 
V,Y,IRES— 1/1/S, N,EPSI $ STATIC SOLUTION 
USET,PG,ULV,UOOV,YS,GO,GM,PS,KFS,KSS,/UGV,PGG,QG/l/ 
*BKL0* S RECOVER DEPENDENT DISPLACEMENTS 
ECT,EPT,BGPDT,SIL,GFfT,CSTM/Xl,X2,X3,ECPT,GPCT/LUSET/ 
NOSIMP/O/NOGENL/GENEL $ TABLES FOR DIFF. STIFF. 
CASECC,GPrr,SIL,EDT,UGV,CSTM,MPT,ECPT,GPCT,DlT/KDGG/ 
S,N,DSCOSET $ DIFF. STIFF. MATRIX 

KDGG,KDNN/MPCF2 / MGG,MNN/MPCF2 $ EQUIV IF NO MPC’S 
LBL1D,MPCF2 $ TRANSFER IF NO MPC’S 
USET,GM,KDGG,,,/KDNN,„ $ MPC’S ON DIFF. STIFF. 

LBL1D S 

KDNN.KDFI7SINGLE / MNN.MFF/SINGLE $ EQUIV. IF NO SPC’S 
FBI, 2D, SINGLE $ TRANSFER IF NO SPC’S 


156 



SCEl 

LABEL 

EQUIV 

COND 

SMP2 

LABEL 

PARAMR 

ADD 

ADD 

EQUIV 

LABEL 

DIAGONAL 

ADD 

ADD 

FRRD 


CHKPNT 
CHKPNT 
EXIT $ 

ENDALTER $ 


USET,KDNN,,,/KDFF,KDFS,KDSS,,, $ SPC’S AND DIFF. STIFF. 
LBL2D $ 

KDFF,KDAA/OMIT / MFF,MAA/OMIT $ EQUIV. IF NO OMITS 
LBL3D,OMIT $ TRANSFER IF NO OMITS 
USET,GO,KDFF/KDAA $ OMITS AND DIFF. STIFF. 

LBL3D $ 

//*SUBC*////MHSPC//HSPC $ NEGATE HYDRO. P 
KDD,KDAA/NEWKDD/(1.0,0.0)/MHSPC $ 
KFS,KDFS/NEWKFS/(1.0,0.0)/MHSPC $ 
NEWKDD,KDD//NEWKFS,KFS $ 

LBL4D $ END OF DIFF. STIFF. EFFECTS (HSP) 
KDD/IDENT/*SQUARE*/0. $ D-SET IDENTITY 
IDENT,/IDM/(1. 0,0.0) $ ANOTHER D-SET IDENTITY 
IDENT,/ZERO/(0.0,0.0) $ D-SET ZERO MATRIX 
CASEXX,USETD,DLT,FRL,GMD,GOD,IDENT,ZERO,IDM,,DIT/ 
UDVF,PSF,PDF,PPF/*DISP*/*DIRECT*/LUSETD/MPCF1/ 
SINGLE/OMIT/NONCUP/FRQSET $ PDF, KDD=I, BDD=0, MDD=i 
MDD,KDD,BDD,PDF,PSF,PPF,EQDYN,USETD,GOD,GMD $ 
KFS,BGPDT,ECT,EQEXIN,GPECT,SIL $ 


The above alter does not depend on whether the fluid is interior or exterior to the structure. 
The Step 2 alters, however, depend on whether an interior fluid is present. For Step 2A, the 
following alter is used: 


ALTER 1 $ NASFIUA STEP 2A, COSMIC 1988 RF8 (REVISED 11/7/89) 

ALTER 2,167 $ REPLACE ALL AFTER ’BEGIN’ AND BEFORE ’END’ 

INPUTT2 /DAT2„„//13 $ INTERNAL FLUID 

INPUTT2 /DAT,„,//11 $ READ SURF MATRIX FROM UT1 

MATPRN DAT,DAT2„, $ 

PARAML DAT//*DMI*/1/S/RIGD $ GET RIGID FLAG 

PARAMR //*EQ*//RIGD/0.////ELAST $ SET ELAST=-1 IF ELASTIC 
COND LBL9D,ELAST $ IF ELASTIC, JUMP OVER RIGID/SOFT 

PARAMR //*EQ*//RIGD/2.////SOFT $ SET SOFT=-l IF SOFT BD. 

COND LBL9E,SOFT $ IF SOFT BOUNDARY, JUMP OVER RIGID 

INPUTT2 /E,PI,VEKC,,//11 $ READ SURF MATRICES FROM UT1 

OUTPUT2 PI,E,„ //-l $ INPUTT2 FILE IS OVER WRITTEN (UT1) 
OUTPUT2 „„ //-9 $ EOF 

CHKPNT DAT,VEKC $ 

EXIT $ 

LABEL LBL9E $ BEGINNING OF SOFT ANALYSIS 

INPUTT2 /CT,PI,VEKC„//11 $ READ SURF MATRICES FROM UT1 

TRNSP CT/C $ 

ADD PI, /MPI/(-l. 0,0.0) $ NEGATE PI 

OUTPUT2 MPI,C„, //-l $ INPUTT2 FILE IS OVER WRITTEN (UT1) 
OUTPUT2 H-9 S EOF 

CHKPNT DAT.VEKC $ 

EXIT $ 

LABEL LBL9D $ BEGINNING OF ELASTIC ANALYSIS 

INPUTT5 /G2,A2,„//14 $ INTERNAL FLUID 
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1NPUTT2 

INPUTT5 

INPUTT2 

MATPRN 

MERGE 

MERGE 

MERGE 

MERGE 

MERGE 

MA'PPRN 

PARAML 

PAR AMR 

PARAMR 

PARAMR 

PARAMR 

PARAMR 

ADD5 

MPYAD 

DECOMP 

FBS 

FBS 

ADD 

ADD 

MPYAD 

MPYAD 

MPYAD 

MPYAD 

MERGE 

MERGE 

MERGE 

EQUIV 

MERGE 

MERGE 

MERGE 

ADD 

EQUIV 

OUTPUT2 

OUTPUT2 

CIIKPNT 

Cl I KPN I 

END ALTER 


/C2T,E2,PI2,VEKC2,//13 $ INTERNAL FLUID 
/G1,A1,,,//12 $ READ SURF MATRICES FROM UT2 
/C1T,E1,PI1,VEKC,FVEC//11 $ READ SURF MATRICES 
EVEC,,,, $ 

A1,,,A2,FVEC,/A/-1 $ 

E1,„E2,FVEC,/E/-1 $ 

C1T,,,C2T,FVEC,/CT/-1 $ 

Gl,,G2,,FVEC,/G/0 $ 

PIl,,,„FVEC/PI/0 $ 

VECM,VECS,,, $ 

DAT//*DMI*/l/2/FREQ $ GET FREQ FROM DA I 

//*complex*//freq/o./freqc $ freq+i*o 

//*MPYC*////W/FRE.QC/(6. 283185,0.) $ OMEGA 
//*MPYC*////IW/W/(0.,1.) $ POMEGA 
//*MPYC*////WW/W/W S OMEGA* *2 
//*SUBC*////MWW//WW S -OMEGA **2 
MDD,KDD,BDD,,/Y/MWW/(1.0,0.0)/IW $ 

G,A,/GA/0 S 

Y/L,U/1//S,N,MINDIAG///S,N,SING $ 

L.U.GA/YIGA/1 $ 

L,U,PDF/YIF/1 $ 

YIGA,/ZIGA/IW $ 

YIF,/ZIF/IW $ 

G,ZIGA,/GTZIGA/1 $ 

CT,GTZIGA,E/II/1 $ LTIS 
G,Z1P,/GTZIF/1 $ 

CT,GTZIF,/Q/'l $ MECHANICAL RHS 
DUM,, PDF,, VECM, /PDF1/1 $ MERGE IN 0 COLUMNS 
DUM , , PSF, , VECM ,/PSFl/'l $ MERGE IN 0 COLUMNS 
DUM, ,PPF, , VECM,/PPF1 /I S MERGE IN 0 COLUMNS 
PDF! , PDF/ /PSF1 , PS F/ /PPF1 , PPF $ 

DUM, ,Q„’ VECM, /RHS1/1 $ MERGE IN ZERO COLUMNS 
DUM, ,GTZIF, , VECM, /GTZIFE/1 S MERGE IN 0 COLUMNS 
DUM,,PI,,VECS,/RIIS2/I $ MERGE IN ZERO COLUMNS 
RHS1,RHS2/RHS $ ADD MECH. AND INC. RFIS 
USETD,DUMl//GOD,DUM2//GMD,DUM3//KFS,DUM4 $ 
RFIS, FI, „ //-l $ INPUIT2 FILE IS OVER WRITTEN (UT1) 

,,,, //-9 $ EOF 

GTZIGA,GTZIFE,GA,PDF,L,U,PSF,DAT,VEKC,FVEC $ 

usf;td,god,gmd,kfs $ 

$ 


The differences between this alter and one used for submerged evacuated structures are due 
to the need to read and combine two sets of SURF matrices, one for each fluid domain. For 
Step 2B, the following alter is used: 


ALTER 1 S NASHUA STEP 2B, COSMIC 1988 RF8 (REVISED 11/7/89) 

ALTER 2,167 $ REPLACE ALL AFTER ’BEGIN’ AND BEFORE ’END’ 

INPUTT2 /PC,,,,// 11 $ READ PRESSURES FROM BLOCK SOLVER (UT1) 
PAR I N PC..I VFC/P1,,,/Q $ REMOVE INTERNAL FLUID DOF 

PAR I N Pl„VEKC/P„,/0 S REMOVE CHIEF DOE FROM P 
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COND 

OUTPUT2 

OUTPUT2 

MATPRN 

EXIT 

LABEL 

MPYAD 

MPYAD 

FBS 

SDR1 

PARTN 

PARTN 

OUTPUT2 

OUTPUT2 

MATPRN 

ENDALTER 


LBL9D,ELAST $ IF ELASTIC, JUMP OVER RIGID/SOFT 
DAT,I>,„ //-l $ INPUTT2 FILE IS OVER-WRITl'EN (UT1) 

„„ //-9 $ EOF 

DA'[',P„, S FOR SOFT BOUNDARY, P REPRESENTS VELOCITY 
$ 

LBL9D $ ELASTIC ANALYSIS 

GTZIGA,PC,GTZIFE/VNC/0/-l $ NORMAL VELOCITIES 
GA,PC,PDF/FA/0/-l $ A-SET FORCES 
L,U,FA/UDVF/1 $ A-SET DISPLACEMENTS 
USETD,,UDVF,,,GOD,GMD,PSF,KFS,,/UPVC,,QPC/l/ 
♦DYNAMICS* $ 

VNC„FVEC/Vl„,/0 $ REMOVE INTERNAL FLUID DOF 
Vl„VEKC/VN,,,/0 $ REMOVE CHIEF DOF FROM VN 
DAT,P,VN,UPVC, U-\ $ INPU1T2 FIFE: IS OVER WRITTEN 
,,„ //- 9 S EOF 
DAT,P,VN,, $ 

$ 


This alter differs from one for evacuated structures because of the presence of several matrix 
partitionings to remove the internal fluid DOF from the solution vectors before the solutions 
are merged with the results for other frequencies. 


NUMERICAL EXAMPLE 

Here we illustrate and validate the two-fluid boundary element formulation developed 
above by solving the problem of acoustic scattering from a submerged fluid-filled spherical 
thin shell. The incident loading is a time-harmonic planar wavctrain, as shown in Fig. 2. 
The specific problem solved has the following characteristics: 


shell mean radius (a) 

5 in 

shell thickness (h) 

0.15 m 

shell Young’s modulus (E) 

2.07 x 10 11 N/m 2 

shell Poisson’s ratio ( i /) 

0.3 

shell density (/> s ) 

7669 kg/m 3 

shell loss factor ( rj ) 

0.01 

fluid density ( p ) 

1000 kg/m 3 

fluid sound speed (c) 

1524 m/s 


The same fluid is used for both the exterior and interior fluid domains. The solution of this 
problem exhibits rotational symmetry about the spherical axis parallel to the direction of 
wave propagation. The benchmark solution to which the numerical results will be compared 
is a scries solution, the derivation of which is summarized in the next section. 


Series Solution 

The series solution for scattering from a submerged evacuated spherical thin shell was 
presented by Junger and Feit (Ref. 13). Here we summarize that solution and indicate the 
modification necessary to include the addition of an interior fluid which fills the spherical 
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Fig. 2. Plane Wave Scattering from a Fluid-Filled Spherical Shell 

volume. 


In general, the series solution for plane wave scattering from a submerged, evacuated, 
spherical thin shell involves computing the impedances of the shell and exterior fluid, the 
scattered field due to rigid body effects, and the radiated field due to elastic (vibrational) 
effects. The shell impedance (the ratio of pressure to normal velocity) for the nth 
axisymmetric shell mode is 


jftc,, h_ [tf-(n< 2 >) 2 ] 

" n a \<f )(„+>, ,-l\ ’ ' * 

where p s is the structural mass density, c p = VE/[p s (l— ir)] , E is Young’s modulus, u is 
Poisson’s ratio, $2= u;a/c p is the dimensionless frequency, h is the shell thickness, a is the 
shell mean radius, /? = h/(a v 12 )^ and X n = n(n+l). The quantities and n( 2) are the 
upper and lower shell resonance dimensionless frequencies, respectively. They are the 
solutions of the characteristic equation 


it -[l+3/y+X n — /9 2 (1— 1/-\ 2 — iA n ] it 

+ (\i — 2)(1— <^+/? 2 [X n 3 — 4X n 2 +\(5— za)— 2(1— tr)] = 0. (19) 


The impedance of the exterior fluid, found by using the Green’s function and identity for the 
exterior fluid, is 


z 


n 


i pc 


h n (ka) 

h'n(ka) ’ 


( 20 ) 


where h n is the Bessel’s function of the third kind of order n. 


Thus, Junger and Feit showed that the far-field scattered pressure is 
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, R>>a, 


( 21 ) 


p(R,*)= 


ie ikR p 0 (2n+l)P n (costf) 
kR ,^0 h' n (ka) 


j'n(ka)- 


££ 

(Z n +z„)(ka) 2 h' n (ka) 


where R is the distance to the field point, 6 is the angle from the z-axis, p c is the incident 
pressure, P n is the Legendre polynomial of order n, and j n is the Bessel’s function of the first 
kind of order n. The two terms in the bracketed expression correspond to rigid body and 
radiated effects, respectively. 


The above expression for the pressure scattered from an evacuated shell can be 
extended to include the effects of the interior fluid merely by replacing the exterior fluid 
impedance z n in Eq. 21 with the sum of the fluid impedances for the exterior and interior 
fluids. It can be shown, by using the Green’s function and identity for the interior domain, 
that the interior impedance, denoted f n , is given by 


Sn 


— i/>c 


jn(ka) 
j'n(ka) ‘ 


( 22 ) 


We note the resemblance between Eqs. 20 and 22 tor the exterior and interior domains, 
respectively. 

The computer program used to evaluate this series solution is a modification of a 
program called SCATSPHERE written by F.M. Henderson, a retired employee of DTRC. 
SCATSPHERE in turn is a variant of an earlier program called RADSPHERE (Ref. 14) for 
computing the radiation from an internally-driven submerged spherical shell. 


Numerical Solution 

A NASTRAN finite element model of the spherical shell was prepared using 40 
axisymmetric conical shell elements spanning the 180 degrees between the two poles of the 
sphere. Due to the axisymmetry of the incident pressure loading, only the N = 0 harmonic 
was required. Since all structural points are in contact with both interior and exterior fluids, 
the resulting model therefore had 205 independent structural degrees of freedom (DOF) and 
41 fluid DOF for each of the two fluid domains. System matrices for the exterior fluid were 
also augmented by the addition of four constraint equations associated with interior Chief 
points to ensure uniqueness of the integral representation at the upper frequencies. The 
nondimensional frequency range 0<ka<5 was swept using a frequency increment of about ka 
= 0.05 with NASHUA and ka = 0.005 with the series solution. Since the series solution is 
converged, we treat it as an "exact” solution for this problem. 

The comparison between the computed and exact solutions is presented is Figs. 3 and 
4, which plot the frequency response of the nondimensional scattered pressure pr/(p 0 a), 
where p is the far-field scattered pressure at distance r from the origin, p Q is the incident 
pressure, and a is the mean radius of the spherical shell. These two figures show very good 
agreement between the two scattering solutions in the backward (0 = 0) and forward (0 - 180 
degrees) directions. In fact, the computed and series solutions are virtually indistinguishable 
from each other. 
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K A 

Fig. 3. Forward Scattering from a Fluid-Filled Spherical Shell 



K A 


Fig. 4. Backward Scattering from a Fluid-Filled Spherical Shell 

DISCUSSION 

A very general computational capability has been described for predicting the sound 
pressure lield radiated or scattered by arbitrary, submerged, fluid-filled, three-dimensional 
elastic structures subjected to time-harmonic loads. The structure is modeled with 
NASTRAN (in all the generality that NASTRAN allows) in combination with boundary 
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element models of both interior and exterior fluid domains. Sufficient automation is 
provided so that, for many structures of practical interest, an existing structural model can 
be adapted for NASHUA acoustic analysis within a few hours. 

One of the many benefits of having NASHUA linked with NASTRAN is the ability to 
integrate the acoustic analysis of a structure with other dynamic analyses. Thus the same 
finite element model can be used for modal analysis, frequency response analysis, linear 
shock analysis, and underwater acoustic analysis. In addition, many of the pre- and 
postprocessors developed for use with NASTRAN become available for NASHUA as well. 
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NASA Ames Research Center Butler Analyses 


In applying a Ritz modal expansion to the solution of a 
transient response, there is a problem as to how many modes are 
needed to obtain accuracy to within a specified percentage. 

One of us, Chargin, has suggested a method based on the 
characteristics of the forcing function. The method can be 
incorporated into the Ritz generation algorithm such that it will 
automatically monitor, regulate and terminate the process 
according to a specified tolerance. 

FORCING FUNCTION CHARACTERISTICS 

Assume that the forcing function F(x,t) can be repre- 
sented as a product of a spatial function and a temporal func- 
tion; i.e. 

(1) F(x, t ) = F(x) . f ( t) . 

Develop a criterion based upon measuring the amount of power 
developed in the forcing function F(x,t). The total power in 
F(x,t) is the product of the "power" in F(x) and f(t) owing to 
the ass um ption in equation (1). The scheme in outline is to 
measure the temporal power and the spatial power in F(x,t) separ- 
ately then compare the corresponding power in the Ritz modes 
against power in each component of the forcing. Generation of 
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additional Ritz inodes will continue until the power criteria are 

met. First a measure is taken of the total power in the temporal 
function. 


Temporal Power 


( 2 ) 


P(x) = 


>t,] 

0 L J 


dt 


where t is the interval over which the transient will act. There 
could very well be a separate temporal function for each spatial 

function. In that case of multiple loadings P(t) of equation (2) 
would be a vector "1" long. 


In order to tailor this power to our use as a guide in 
selecting Ritz modes it will be useful to measure the amount of 
temporal power as a function of frequency. Expand the temporal 

function in a Fourier Series and sum the power versus the expan- 
sion multiple: 


(3) f(t) = a Q + a^cos wt + b-^sin wt + a 2 cos 2wt + b 2 sin 2wt 

+ +a cos wt + b sin wt . 

n n 

The power within a band 0 to nw is: 

_ n ( . 2 _n 7 

(4) P = > a. cos iwt + b.sin iwt = > af + bf . 

n 1 = 0 ^ 1 1 ' 1=0 1 x > 

One can compare the amount of power within a given band P with 
the total power P(x). When P^ is within close range of P(x), say 
1%, then the analyst can be satisfied that the frequency range of 
the truncated temporal forcing function is sufficiently broad to 
encompass the temporal requirements of the forcing. 


166 


MONITORING OF RIT2 MODAL GENERATION 



p 

For a certain frequency f n = ^ = ^et this fre- 

quency be f . 

Spatial Power 

Now turn to the spatial distribution of the forcinq function 
F(x), and develop a measure that is called spatial power. Make 
an additional assumption that all points being loaded have mass. 

(5) CTI(x)3 = CF(x)3 T CMD CF(x) 3 = [c ia ] [f^>] ' 

where [ c i a ] = |V(x) T mJ is a coefficient matrix that will be shown 

to be usefual later. II(x) is a matrix if there are a number of 
loading cases "1". 
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At this point we have 2 measures of the forcing fuction. 
We have f^ reflecting the desired frequency content and IT(x) 
reflecting the desired level of power to activate the structural 
mass. Now it is time to test the adequacy of the number of Ritz 
modes to repond to the forcing at these power levels. 

First, get a spatial measure of the Ritz modes 0^. A 
simple scheme is to pattern the measure after equation (5) but 
substitute the matrix of k spatial Ritz functions for the post- 
multiply operation instead of the spatial component of loading 
F( x) . 


<6, ['R rv ] = ftF«x.j T c M] [e k ] 

The way to use R^ is to compare each diagonal term of H(x) with 
the corresponding diagonal term of to :ee if the ratio is 
within a specified tolerance; i.e. 


(7) 



< e for every 1. 


Keep generating additional Ritz modes 0^ k > r until an r has 
been reached for which every diagonal term satisfies the 
criterion . 


Monitoring Function 

If the 0^ satisfy the spatial power criterion of F(x,t) it does 
not necessarily hold that 0^ will simultaneously satisfy the 
temporal requirements of F(x,t). Therefore, use the 0^ which 
have met the spatial power requirements and obtain an estimate of 
its frequency content by setting up the frequency equation. Use 
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an estimator that is much less demanding' than solving '-he 
eigenvalue problem. Merely do the following. 

Compute a k— order Ritz generalized mass and stiffness; i.e. 


( 8 ) 



and kJ = |_ e k_ 


and construct a test 
threshold frequency, f , 


matrix called 
determined from 


[S k ] which involves the 
the temporal power P . 


(9) [k^ - 5 [Sj. 

Decompose and extract the value of the output parameter 

NBRCHG issued by the DECOMP module. Param NBRCHG reports the 

number of negative values on the factor diagonal of which is 

tantamount to the number of sign changes or zero crossings in the 

characteristic equation. It the value of NBRCHG = k, it 

indicates that the frequencies of all k Ritz modes are less than 

f . One would be inclined to want the frequency content of the 
o 

Ritz modes to bracket f ; i.e. some modes with a frequecy ; t Q . 
This implies that one would seek to have the value of NERCHG to 
be less than the order of matrix S^. This idea can be built into 
the Ritz generation routine as a test as to whether enough modes 
have been generated to within a certain margin ot such that 
af Q , where the user specifies a. 


Obviously, one would not want to repeat the eigenvalue 

T 

estimate each time a new Ritz mode is obtained, because 
and B^M 0^ could become expensive as k becomes large. Therefore 
develop a scheme whereby the spatial power is reduced by some 
factor y > 1 ; i.e. e = e , . / y • Typically, one can use a value 

1 X C W w -L V— ^ 


169 


MONITORING OF RITZ MODAL GENERATION 


of 2 to 4 for y. Once the spatial error e is satisried, 

new 

repeat the eigenvalue estimate test. 

Conclusion 


A scheme has been proposed to monitor the adequacy of a 
set of Ritz modes to represent a solution by comparing the 
quantity generated with certain properties involving the forcing 
function. In so doing an attempt has been made to keep this 
algorithm lean and efficient, so that it will be economical to 
apply. Using this monitoring scheme during Ritz Mode generation 
will automatically ensure that the k Ritz modes 9^ that are 
generated are adequate to represent both the spatial and temporal 
behavior of the structure when forced under the given transient 
condition defined by F(x,t). 
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